#' Estimate Person Parameters
#'
#' Compute person parameters for the 1,2,3,4-PL model and for the GPCM. Choose between ML, WL, MAP, EAP and robust estimation. Use this function if 4-PL items and GPCM items are mixed for each person.
#'
#'
#' For a test with both: dichotomous and polytomous items which have been scaled under 1/2/3/4-PL model or the GPCM, use this function to estimate the person ability parameters. You have to define the appropriate model for each item.
#'
#'Please note, that \code{robust} estimation with (Huber ability estimate) polytomous items is still experimental!
#'
#'@param respm An integer matrix, which contains the examinees responses. A persons x items matrix is expected.
#'@param thres A numeric matrix which contains the threshold parameter for each item. If the first row of the matrix is not set to zero (only zeroes in the first row) - then a row-vector with zeroes is added by default.
#'@param slopes A numeric vector, which contains the slope parameters for each item - one parameter per item is expected.
#'@param lowerA A numeric vector, which contains the lower asymptote parameters (kind of guessing parameter) for each item. In the case of polytomous items, the value must be 0.
#'@param upperA numeric vector, which contains the upper asymptote parameters for each item. In the case of polytomous items, the value must be 1.
#'@param theta_start A vector which contains a starting value for each person. If NULL is submitted, the starting values are set automatically. If a scalar is submitted, this start value is used for each person.
#'@param mu A numeric vector of location parameters for each person in case of MAP estimation. If nothing is submitted this is set to 0 for each person for MAP estimation.
#'@param sigma2 A numeric vector of variance parameters for each person in case of MAP or EAP estimation. If nothing is submitted this is set to 1 for each person for MAP estimation.
#'@param type Which maximization should be applied? There are five valid entries possible: "mle", "wle", "map", "eap" and "robust". To choose between the methods, or just to get a deeper understanding the papers mentioned below are quite helpful. The default is \code{"wle"} which is a good choice in many cases.
#'
#'@param model2est A character vector with length equal to the number of submitted items. It defines itemwise the response model under which the item parameter was estimated. There are 2 valid inputs up to now: \code{"GPCM"} and \code{"4PL"}.
#'
#'@param maxsteps The maximum number of steps the NR algorithm will take. Default = 100.
#'@param exac How accurate are the estimates supposed to be? Default is 0.001.
#'@param H In case \code{type = "robust"} a Huber ability estimate is performed, and \code{H} modulates how fast the downweighting takes place (for more Details read Schuster & Yuan 2011).
#'@param ctrl More controls:
#'
#'\itemize{
#' \item \code{killdupli}: Should duplicated response pattern be removed for estimation (estimation is faster)? This is especially resonable in case of a large number of examinees and a small number of items. Use this option with caution (for map and eap), because persons with different \code{mu} and \code{sigma2} will have different ability estimates despite they responded identically. Default value is \code{FALSE}.
#'
#'\item \code{skipcheck}: Default = FALSE. If TRUE data matrix and arguments are not checked - this saves time e.g. when you use this function for simulations.
#'
#'}
#'
#'@template resulttemplate
#'
#' @seealso \link{PPass}, \link{PP_gpcm}, \link{PP_4pl}, \link{JKpp}, \link{PV}
#'
#'@export
#'
#'@author Manuel Reif
#'
#'@references Baker, Frank B., and Kim, Seock-Ho (2004). Item Response Theory - Parameter Estimation Techniques. CRC-Press.
#'
#'Barton, M. A., & Lord, F. M. (1981). An Upper Asymptote for the Three-Parameter Logistic Item-Response Model.
#'
#'Magis, D. (2013). A note on the item information function of the four-parameter logistic model. Applied Psychological Measurement, 37(4), 304-315.
#'
#'Muraki, Eiji (1992). A Generalized Partial Credit Model: Application of an EM Algorithm. Applied Psychological Measurement, 16, 159-176.
#'
#'Muraki, Eiji (1993). Information Functions of the Generalized Partial Credit Model. Applied Psychological Measurement, 17, 351-363.
#'
#'Samejima, Fumiko (1993). The bias function of the maximum likelihood estimate of ability for the dichotomous response level. Psychometrika, 58, 195-209.
#'
#'Samejima, Fumiko (1993). An approximation of the bias function of the maximum likelihood estimate of a latent variable for the general case where the item responses are discrete. Psychometrika, 58, 119-138.
#'
#'Schuster, C., & Yuan, K. H. (2011). Robust estimation of latent ability in item response models. Journal of Educational and Behavioral Statistics, 36(6), 720-735.
#'
#'Wang, S. and Wang, T. (2001). Precision of Warm's Weighted Likelihood Estimates for a Polytomous Model in Computerized Adaptive Testing. Applied Psychological Measurement, 25, 317-331.
#'
#'Warm, Thomas A. (1989). Weighted Likelihood Estimation Of Ability In Item Response Theory. Psychometrika, 54, 427-450.
#'
#'Yen, Y.-C., Ho, R.-G., Liao, W.-W., Chen, L.-J., & Kuo, C.-C. (2012). An empirical evaluation of the slip correction in the four parameter logistic models with computerized adaptive testing. Applied Psychological Measurement, 36, 75-87.
#'
#'@example ./R/.examples.R
#'@rdname PPall
#'
PPall <- function(respm, thres, slopes, lowerA, upperA, theta_start=NULL,
mu = NULL, sigma2 = NULL, type="wle", model2est,
maxsteps=100, exac=0.001, H=1, ctrl=list())
{
###
call <- match.call()
attr(call, "date") <- date()
attr(call,"version") <- packageVersion("PP")
###
## --------- user controls
cont <- list(killdupli=FALSE,skipcheck=FALSE)
user_ctrlI <- match(names(ctrl),names(cont))
if(any(is.na(user_ctrlI)))
{
notex <- names(ctrl)[is.na(user_ctrlI)]
warning("The following options in ctrl do not exist: ", paste(notex,collapse=", "))
ctrl <- ctrl[!is.na(user_ctrlI)]
user_ctrlI <- user_ctrlI[!is.na(user_ctrlI)]
}
cont[user_ctrlI] <- ctrl
######## take care of the threshold matrix
if(is.matrix(thres))
{
#iimm <- nrow(thres) == 1
if(any(thres[1,] != 0))
{
thres <- rbind(0,thres)
}
# compute the maximal score per item
maxsc <- apply(thres,2,function(x)(length(x) - sum(is.na(x)))-1)
# are there any items with more than 2 categories?
#allebigger <- any(maxsc > 1)
}
## starting values
if(is.null(theta_start))
{
thresv <- apply(thres[2,,drop=FALSE],2, function(tr) mean(tr,na.rm=TRUE))
theta_start <- apply(respm, 1, function(st)
{
mean(thresv[!is.na(st)])
})
} else if(length(theta_start) == 1)
{
theta_start <- rep(theta_start, nrow(respm))
}
## ----------------------
# check for errors and warnings etc --------------
## --------- check user inputs
if(!cont$skipcheck) # to save time e.g. when running simulations
{
## --------- check user inputs
checkINP(respm, thres, slopes, theta_start, lowerA, upperA, type)
# -----------------------------------------------------
}
model2est <- match.arg(model2est,c("GPCM","4PL"),several.ok = TRUE)
if(length(model2est) < ncol(respm)) stop("Check the properties for the submitted model2est vector!\n")
# in case map is chosen and no mu and/or sigma2 is/are submitted.
if( (any(is.null(mu)) | any(is.null(sigma2))))
{
if(any(is.null(mu)))
{
mu <- rep(0,nrow(respm))
if(type == "map") warning("all mu's are set to 0! \n")
}
if(any(is.null(sigma2)))
{
sigma2 <- rep(1,nrow(respm))
if(type == "map") warning("all sigma2's are set to 1! \n")
}
}
## ----------------------
# ----- MLE: Inf and NA ---------------#
if(type=="mle" | type=="robust")
{
resPPx <- ansol(respm,maxsc)
respm <- respm[!is.na(resPPx[,2]),]
}
# ----- remove duplicated -------------#
if(cont$killdupli)
{
# this vector contains the indices, to recreate the large matrix in the end of the estimation
dupvec <- make_dup(respm)
respm <- respm[dupvec$ndpat,,drop=FALSE]
}
# ----- estimation procedure -------------#
if(any(lowerA[model2est == "GPCM"] != 0)) warning("At least for one item with #categories > 1, there is an lower asymptote != 0 submitted. This will be ignored.\n")
if(any(upperA[model2est == "GPCM"] != 1)) warning("At least for one item with #categories > 1, there is an upper asymptote != 1 submitted. This will be ignored.\n")
modelcat <- paste(unique(model2est),collapse=", ")
cat("Estimating: mixed ",modelcat, "... \n")
cat("type =",type,"\n")
if(type %in% c("mle","wle","map"))
{
# ----- estimation procedure -------------#
resPP <- NR_mixed(respm,DELTA = thres,ALPHA = slopes, LOWA = lowerA, UPPA = upperA, THETA = theta_start,model=model2est, wm=type,maxsteps,exac,mu,sigma2,H)
resPP$resPP[,2] <- sqrt(resPP$resPP[,2])
} else if(type == "robust")
{
warning("Robust estimation for GPCM is still very experimental! \n")
resPP <- NR_mixed(respm,DELTA = thres,ALPHA = slopes, LOWA= lowerA, UPPA = upperA, THETA = theta_start,model=model2est, wm=type,maxsteps,exac,mu,sigma2,H=H)
resPP$resPP[,2] <- sqrt(resPP$resPP[,2])
} else if(type == "eap")
{
resPP <- list()
resPP$resPP <- eap_mixed(respm, thres, slopes, lowerA=lowerA,
upperA=upperA, mu = mu, sigma2 = sigma2, model2est=model2est)
resPP$nsteps <- 0
}
### result preperation --------------------------
if(cont$killdupli)
{
resPP$resPP <- resPP$resPP[dupvec$posvec,]
}
if(type=="mle" | type=="robust")
{
resPPx[!is.na(resPPx[,2]),] <- resPP$resPP
resPP$resPP <- resPPx
}
## ---------------------------------------------
colnames(resPP$resPP) <- c("estimate","SE")
ipar <- list(respm=respm,thres=thres,slopes=slopes,lowerA=lowerA,
upperA=upperA,theta_start=theta_start,mu=mu,sigma2=sigma2,
cont=cont,model2est=model2est,H=H)
if(cont$killdupli)
{
ipar$dupvec <- dupvec
}
## ---------------------------------------------
cat("Estimation finished!\n")
rescall <- list(resPP=resPP,call=call,type=type,ipar=ipar)
class(rescall) <- c("gpcm4pl","ppeo")
return(rescall)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.