Nothing
R/SingleGroup-methods.R
In mirt: Multidimensional Item Response Theory
Defines functions
traditional2mirt
mirt2traditional
Documented in
traditional2mirt
#' Print the model objects
#'
#' Print model object summaries to the console.
#'
#' @param x an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{MixedClass}
#'
#' @name print-method
#' @aliases print,SingleGroupClass-method
#' print,MultipleGroupClass-method print,MixedClass-method print,DiscreteClass-method
#' print,MixtureClass-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @rdname print-method
#' @examples
#'
#' \dontrun{
#' x <- mirt(Science, 1)
#' print(x)
#' }
setMethod(
f = "print",
signature = signature(x = 'SingleGroupClass'),
definition = function(x){
if(!length(x@time)){
cat('An object of class \"SingleGroupClass\"\n')
return(invisible(NULL))
}
cat("\nCall:\n", paste(deparse(x@Call), sep = "\n", collapse = "\n"),
"\n\n", sep = "")
cat("Full-information item factor analysis with ", x@Model$nfact, " factor(s).\n", sep="")
EMquad <- ''
if(x@Options$method == 'EM') EMquad <- c('\n using ', x@Options$quadpts, ' quadrature')
method <- x@Options$method
if(method == 'MIXED') method <- 'MHRM'
if(x@OptimInfo$converged)
cat("Converged within ", x@Options$TOL, ' tolerance after ', x@OptimInfo$iter, ' ',
method, " iterations.\n", sep = "")
else
cat("FAILED TO CONVERGE within ", x@Options$TOL, ' tolerance after ',
x@OptimInfo$iter, ' ', method, " iterations.\n", sep="")
cat('mirt version:', as.character(utils::packageVersion('mirt')), '\n')
cat('M-step optimizer:', x@Options$Moptim, '\n')
if(method %in% c('EM', 'QMCEM', 'BL', 'MCEM')){
if(method == 'EM' || method == 'QMCEM')
cat('EM acceleration:', x@Options$accelerate, '\n')
if(method == 'EM' || method == 'BL')
cat('Number of rectangular quadrature:', x@Options$quadpts)
else if(method == 'QMCEM')
cat('Number of quasi-Monte Carlo points:', x@Options$quadpts)
else if(method == 'MCEM')
cat('Number of Monte Carlo points:', x@Options$quadpts)
cat('\n')
}
dentype <- switch(x@Options$dentype,
"EH" = "Empirical histogram",
"EHW" = 'Empirical histogram (scaled)',
x@Options$dentype)
cat('Latent density type:', dentype, '\n')
if(method == 'MHRM')
cat("Average MH acceptance ratio(s):", paste0(round(x@OptimInfo$aveAR,3), collapse=', '), "\n")
if(!is.na(x@OptimInfo$secondordertest)){
cat("\nInformation matrix estimated with method:", x@Options$SE.type)
cat('\nSecond-order test: model ', if(!x@OptimInfo$secondordertest)
'is not a maximum or the information matrix is too inaccurate' else
'is a possible local maximum', sep = "")
if(x@OptimInfo$secondordertest)
cat("\nCondition number of information matrix = ", x@OptimInfo$condnum)
cat('\n')
}
if(length(x@Fit$logLik) > 0){
if(x@Fit$logPrior != 0){
cat("\nLog-posterior = ", x@Fit$logLik + x@Fit$logPrior, if(method == 'MHRM')
paste(', SE =', round(x@Fit$SElogLik,3)), "\n",sep='')
cat('Estimated parameters:', length(extract.mirt(x, 'parvec')), '\n')
} else {
cat("\nLog-likelihood = ", x@Fit$logLik, if(method == 'MHRM')
paste(', SE =', round(x@Fit$SElogLik,3)), "\n",sep='')
cat('Estimated parameters:', extract.mirt(x, 'nestpars'), '\n')
}
cat("AIC = ", x@Fit$AIC, "\n", sep='')
cat("BIC = ", x@Fit$BIC, "; SABIC = ", x@Fit$SABIC, "\n", sep='')
if(!is.nan(x@Fit$p)){
cat("G2 (", x@Fit$df,") = ", round(x@Fit$G2,2), ", p = ", round(x@Fit$p,4), sep='')
cat("\nRMSEA = ", round(x@Fit$RMSEA,3), ", CFI = ", round(x@Fit$CFI,3),
", TLI = ", round(x@Fit$TLI,3), sep='')
}
}
cat("\n")
}
)
#' Show model object
#'
#' Print model object summaries to the console.
#'
#' @param object an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{MixedClass}
#'
#' @name show-method
#' @aliases show,SingleGroupClass-method
#' show,MultipleGroupClass-method show,MixedClass-method show,DiscreteClass-method
#' show,MixtureClass-method
#' @rdname show-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @examples
#'
#' \dontrun{
#' x <- mirt(Science, 1)
#' show(x)
#' }
setMethod(
f = "show",
signature = signature(object = 'SingleGroupClass'),
definition = function(object){
print(object)
}
)
#' Summary of model object
#'
#' Transforms coefficients into a standardized factor loading's metric. For \code{MixedClass} objects,
#' the fixed and random coefficients are printed. Note that while the output to the console is rounded
#' to three digits, the returned list of objects is not. For simulations, use
#' \code{output <- summary(mod, verbose = FALSE)} to suppress the console messages.
#'
#' @param object an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{MixedClass}
#' @param rotate a string indicating which rotation to use for exploratory models, primarily
#' from the \code{GPArotation} package (see documentation therein).
#'
#' Rotations currently supported are: \code{'promax'}, \code{'oblimin'}, \code{'varimax'},
#' \code{'quartimin'}, \code{'targetT'}, \code{'targetQ'}, \code{'pstT'}, \code{'pstQ'},
#' \code{'oblimax'}, \code{'entropy'}, \code{'quartimax'}, \code{'simplimax'},
#' \code{'bentlerT'}, \code{'bentlerQ'}, \code{'tandemI'}, \code{'tandemII'},
#' \code{'geominT'}, \code{'geominQ'}, \code{'cfT'}, \code{'cfQ'}, \code{'infomaxT'},
#' \code{'infomaxQ'}, \code{'mccammon'}, \code{'bifactorT'}, \code{'bifactorQ'}.
#'
#' For models that are not exploratory this input will automatically be set to \code{'none'}
#' @param Target a dummy variable matrix indicting a target rotation pattern. This is required for
#' rotations such as \code{'targetT'}, \code{'targetQ'}, \code{'pstT'}, and \code{'pstQ'}
#' @param suppress a numeric value indicating which (possibly rotated) factor
#' loadings should be suppressed. Typical values are around .3 in most
#' statistical software. Default is 0 for no suppression
#' @param suppress.cor same as \code{suppress}, but for the correlation matrix
#' output
#' @param verbose logical; allow information to be printed to the console?
#' @param ... additional arguments to be passed
#'
#' @name summary-method
#' @aliases summary,SingleGroupClass-method
#' summary,MultipleGroupClass-method summary,MixedClass-method summary,DiscreteClass-method
#' summary,MixtureClass-method
#' @export
#' @rdname summary-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @seealso \code{\link{coef-method}}
#' @examples
#'
#' \dontrun{
#' x <- mirt(Science, 2)
#' summary(x)
#' summary(x, rotate = 'varimax')
#'
#' }
setMethod(
f = "summary",
signature = 'SingleGroupClass',
definition = function(object, rotate = 'oblimin', Target = NULL,
suppress = 0, suppress.cor = 0,
verbose = TRUE, ...){
if (!object@Options$exploratory || rotate == 'none') {
F <- object@Fit$F
F[abs(F) < suppress] <- NA
h2 <- as.matrix(object@Fit$h2)
SS <- apply(F^2,2,sum)
gp <- ExtractGroupPars(object@ParObjects$pars[[object@Data$nitems + 1L]])
Phi <- cov2cor(gp$gcov)
colnames(h2) <- "h2"
rownames(Phi) <- colnames(Phi) <- names(SS) <-
colnames(F)[seq_len(object@Model$nfact)]
loads <- cbind(F,h2)
if(verbose){
if(object@Options$exploratory)
cat("\nUnrotated factor loadings: \n\n")
print(loads, 3, na.print = " ")
cat("\nSS loadings: ", round(SS, 3), "\n")
cat("Proportion Var: ",round(SS/nrow(F), 3), "\n")
cat("\nFactor correlations: \n\n")
Phiprint <- Phi
Phiprint[abs(Phi) < suppress.cor] <- NA
Phiprint[upper.tri(Phiprint, diag = FALSE)] <- NA
print(round(Phiprint, 3), na.print = " ")
}
invisible(list(rotF=F,h2=h2,fcor=Phi))
} else {
F <- object@Fit$F
h2 <- as.matrix(object@Fit$h2)
colnames(h2) <- "h2"
rotF <- Rotate(F, rotate, Target = Target, ...)
SS <- apply(rotF$loadings^2,2,sum)
L <- rotF$loadings
L[abs(L) < suppress] <- NA
loads <- cbind(L,h2)
Phi <- diag(ncol(F))
if(!rotF$orthogonal){
Phi <- rotF$Phi
}
colnames(Phi) <- rownames(Phi) <- colnames(F)
if(verbose){
cat("\nRotation: ", rotate, "\n")
cat("\nRotated factor loadings: \n\n")
print(loads, 3, na.print = " ")
cat("\nRotated SS loadings: ",round(SS,3), "\n")
cat("\nFactor correlations: \n\n")
Phiprint <- Phi
Phiprint[abs(Phi) < suppress.cor] <- NA
Phiprint[upper.tri(Phiprint, diag = FALSE)] <- NA
print(round(Phiprint, 3), na.print = " ")
}
if(any(h2 > 1))
warning("Solution has Heywood cases. Interpret with caution.",
call.=FALSE)
invisible(list(rotF=rotF$loadings,h2=h2,fcor=Phi))
}
}
)
#' Extract raw coefs from model object
#'
#' Return a list (or data.frame) of raw item and group level coefficients. Note that while
#' the output to the console is rounded to three digits, the returned list of objects is not.
#' Hence, elements from \code{cfs <- coef(mod); cfs[[1]]} will contain the non-rounded results (useful
#' for simulations).
#'
#' @param object an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{MixedClass}
#' @param CI the amount of converged used to compute confidence intervals; default is
#' 95 percent confidence intervals
#' @param IRTpars logical; convert slope intercept parameters into traditional IRT parameters?
#' Only applicable to unidimensional models or models with simple structure (i.e., only one non-zero slope).
#' If a suitable ACOV estimate was computed in the fitted
#' model, and \code{printSE = FALSE}, then suitable CIs will be included based on the delta
#' method (where applicable)
#' @param rotate see \code{summary} method for details. The default rotation is \code{'none'}
#' @param Target a dummy variable matrix indicting a target rotation pattern
#' @param printSE logical; print the standard errors instead of the confidence intervals? When
#' \code{IRTpars = TRUE} then the delta method will be used to compute the associated standard errors
#' from mirt's default slope-intercept form
#' @param as.data.frame logical; convert list output to a data.frame instead?
#' @param simplify logical; if all items have the same parameter names (indicating they are
#' of the same class) then they are collapsed to a matrix, and a list of length 2 is returned
#' containing a matrix of item parameters and group-level estimates
#' @param unique return the vector of uniquely estimated parameters
#' @param verbose logical; allow information to be printed to the console?
#' @param rawug logical; return the untransformed internal g and u parameters?
#' If \code{FALSE}, g and u's are converted with the original format along with delta standard errors
#' @param ... additional arguments to be passed
#'
#' @name coef-method
#' @aliases coef,SingleGroupClass-method
#' coef,MultipleGroupClass-method coef,MixedClass-method coef,DiscreteClass-method
#' coef,MixtureClass-method
#' @rdname coef-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @export
#' @seealso \code{\link{summary-method}}
#' @examples
#'
#' \dontrun{
#' dat <- expand.table(LSAT7)
#' x <- mirt(dat, 1)
#' coef(x)
#' coef(x, IRTpars = TRUE)
#' coef(x, simplify = TRUE)
#'
#' #with computed information matrix
#' x <- mirt(dat, 1, SE = TRUE)
#' coef(x)
#' coef(x, printSE = TRUE)
#' coef(x, as.data.frame = TRUE)
#'
#' #two factors
#' x2 <- mirt(Science, 2)
#' coef(x2)
#' coef(x2, rotate = 'varimax')
#'
#' }
setMethod(
f = "coef",
signature = 'SingleGroupClass',
definition = function(object, CI = .95, printSE = FALSE, rotate = 'none', Target = NULL,
IRTpars = FALSE, rawug = FALSE, as.data.frame = FALSE,
simplify = FALSE, unique = FALSE, verbose = TRUE, ...){
dots <- list(...)
discrete <- ifelse(is.null(dots$discrete), FALSE, TRUE)
if(unique) return(extract.mirt(object, 'parvec'))
if(printSE && length(object@ParObjects$pars[[1L]]@SEpar)) rawug <- TRUE
if(CI >= 1 || CI <= 0)
stop('CI must be between 0 and 1', call.=FALSE)
z <- abs(qnorm((1 - CI)/2))
SEnames <- paste0('CI_', c((1 - CI)/2*100, ((1 - CI)/2 + CI)*100))
J <- object@Data$nitems
nfact <- object@Model$nfact + length(object@Model$prodlist)
a <- matrix(0, J, nfact)
for(i in 1:J)
a[i, ] <- ExtractLambdas(object@ParObjects$pars[[i]])
if (object@Options$exploratory && rotate != 'none'){
if(verbose) cat("\nRotation: ", rotate, "\n\n")
so <- summary(object, rotate=rotate, Target=Target, verbose=FALSE, ...)
a <- rotateLambdas(so) * 1.702
for(i in 1:J){
object@ParObjects$pars[[i]]@par[1:nfact] <- a[i, ]
object@ParObjects$pars[[i]]@SEpar <- numeric(0L)
}
object@ParObjects$pars[[J + 1]]@par[-c(1:nfact)] <- so$fcor[lower.tri(so$fcor, TRUE)]
}
if(IRTpars){
if(object@Model$nfact > 1L){
apars <- lapply(object@ParObjects$pars[-(extract.mirt(object, 'nitems') + 1L)],
function(x) x@par[1L:object@Model$nfact])
is_ss <- sapply(apars, function(x) sum(x != 0) == 1L)
if(!any(is_ss)){
warning(c('Traditional parameterization only available for unidimensional ',
'models or models with simple structure patterns (neither found)'), call.=FALSE)
IRTpars <- FALSE
}
}
if(IRTpars){
vcov <- vcov(object)
for(i in 1L:J){
if(class(object@ParObjects$pars[[i]]) %in% c('gpcmIRT')) next
object@ParObjects$pars[[i]] <- mirt2traditional(object@ParObjects$pars[[i]],
vcov=vcov, nfact=object@Model$nfact)
}
}
}
allPars <- list()
if(length(object@ParObjects$pars[[1L]]@SEpar) && !simplify){
if(printSE){
for(i in seq_len(J+1L)){
allPars[[i]] <- matrix(c(object@ParObjects$pars[[i]]@par,
object@ParObjects$pars[[i]]@SEpar),
2, byrow = TRUE)
rownames(allPars[[i]]) <- c('par', 'SE')
nms <- names(object@ParObjects$pars[[i]]@est)
if(i <= J && object@Model$itemtype[i] != 'custom' && !IRTpars){
nms[nms == 'g'] <- 'logit(g)'
nms[nms == 'u'] <- 'logit(u)'
}
colnames(allPars[[i]]) <- nms
}
} else {
for(i in seq_len(J+1L)){
allPars[[i]] <- matrix(c(object@ParObjects$pars[[i]]@par,
object@ParObjects$pars[[i]]@par -
z*object@ParObjects$pars[[i]]@SEpar,
object@ParObjects$pars[[i]]@par +
z*object@ParObjects$pars[[i]]@SEpar),
3, byrow = TRUE)
rownames(allPars[[i]]) <- c('par', SEnames)
colnames(allPars[[i]]) <- object@ParObjects$pars[[i]]@parnames
}
}
} else {
for(i in seq_len(J+1L)){
allPars[[i]] <- matrix(object@ParObjects$pars[[i]]@par, 1L)
colnames(allPars[[i]]) <- object@ParObjects$pars[[i]]@parnames
rownames(allPars[[i]]) <- 'par'
}
}
if(!rawug && !IRTpars){
allPars <- lapply(allPars, function(x, digits){
x[ , colnames(x) %in% c('g', 'u')] <-
antilogit(x[ , colnames(x) %in% c('g', 'u')])
x
})
}
names(allPars) <- c(colnames(object@Data$data), 'GroupPars')
if(as.data.frame)
allPars <- t(as.data.frame(allPars))
if(simplify && !as.data.frame){
allPars <- lapply(allPars, function(x) x[1L, , drop=FALSE])
items.old <- allPars[seq_len(length(allPars)-1L)]
nms <- lapply(items.old, colnames)
unms <- unique(do.call(c, nms))
items <- matrix(NA, length(items.old), length(unms))
rownames(items) <- names(items.old)
colnames(items) <- unms
for(i in seq_len(nrow(items)))
items[i, nms[[i]]] <- items.old[[i]]
nfact <- object@Model$nfact
if(discrete){
allPars <- list(items=items, group.intercepts=allPars$GroupPars)
} else if(object@ParObjects$pars[[J+1L]]@dentype != "custom"){
means <- allPars$GroupPars[seq_len(nfact)]
if(object@ParObjects$pars[[J+1L]]@dentype == "Davidian"){
covs <- matrix(NA, nfact, nfact)
covs[lower.tri(covs, TRUE)] <- allPars$GroupPars[2L]
covs[upper.tri(covs, FALSE)] <- covs[lower.tri(covs, FALSE)]
colnames(covs) <- rownames(covs) <- names(means) <-
object@Model$factorNames[seq_len(nfact)]
allPars <- list(items=items, means=means, cov=covs,
Davidian_phis=allPars$GroupPars[-c(1:2)])
} else {
covs <- matrix(NA, nfact, nfact)
if(object@ParObjects$pars[[J+1L]]@dentype == "mixture")
covs[lower.tri(covs, TRUE)] <-
allPars$GroupPars[-c(seq_len(nfact), length(allPars$GroupPars))]
else covs[lower.tri(covs, TRUE)] <- allPars$GroupPars[-seq_len(nfact)]
covs <- makeSymMat(covs)
colnames(covs) <- rownames(covs) <- names(means) <-
object@Model$factorNames[seq_len(nfact)]
allPars <- list(items=items, means=means, cov=covs)
}
} else allPars <- list(items=items, GroupPars=allPars$GroupPars)
}
if(.hasSlot(object@Model$lrPars, 'beta')){
allPars$lr.betas <- object@Model$lrPars@beta
if(!all(is.nan(object@Model$lrPars@SEpar))){
tmp <- allPars$lr.betas
if(printSE){
tmp[] <- object@Model$lrPars@SEpar
allPars$lr.betas <- list(betas=allPars$lr.betas, SE=tmp)
} else {
low <- tmp - z*object@Model$lrPars@SEpar
high <- tmp + z*object@Model$lrPars@SEpar
allPars$lr.betas <- list(betas=allPars$lr.betas, low, high)
names(allPars$lr.betas) <- c('betas', SEnames)
}
}
}
if(!as.data.frame)
class(allPars) <- c('mirt_list', 'list')
return(allPars)
}
)
#' Compare nested models with likelihood-based statistics
#'
#' Compare nested models using likelihood ratio test (X2), Akaike Information Criterion (AIC),
#' Bayesian Information Criterion (BIC),
#' Sample-Size Adjusted BIC (SABIC), and Hannan-Quinn (HQ) Criterion.
#' When given a sequence of objects, \code{anova} tests the models against one another
#' in the order specified. Note that the \code{object} inputs should be ordered in terms
#' of most constrained model to least constrained.
#'
#' @param object an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{MixedClass}, reflecting the most constrained model fitted
#' @param object2 a second model estimated from any of the mirt package estimation methods
#' @param ... additional less constrained model objects to be compared
#' sequentially to the previous model
#' @param bounded logical; are the two models comparing a bounded parameter (e.g., comparing a single
#' 2PL and 3PL model with 1 df)? If \code{TRUE} then a 50:50 mix of chi-squared distributions
#' is used to obtain the p-value
#' @param mix proportion of chi-squared mixtures. Default is 0.5
#' @param frame (internal parameter not for standard use)
#' @param verbose (deprecated argument)
#'
#' @return a \code{data.frame}/\code{mirt_df} object
#'
#' @name anova-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @aliases anova,SingleGroupClass-method
#' anova,MultipleGroupClass-method anova,MixedClass-method anova,DiscreteClass-method
#' anova,MixtureClass-method
#' @export
#' @rdname anova-method
#' @examples
#'
#' \dontrun{
#' x <- mirt(Science, 1)
#' x2 <- mirt(Science, 2)
#' anova(x, x2)
#'
#' # compare three models sequentially (X2 not always meaningful)
#' x3 <- mirt(Science, 1, 'gpcm')
#' x4 <- mirt(Science, 1, 'nominal')
#' anova(x, x2, x3, x4)
#'
#' # in isolation
#' anova(x)
#'
#' # with priors on first model
#' model <- "Theta = 1-4
#' PRIOR = (1-4, a1, lnorm, 0, 10)"
#' xp <- mirt(Science, model)
#' anova(xp, x2)
#' anova(xp)
#'
#' # bounded parameter
#' dat <- expand.table(LSAT7)
#' mod <- mirt(dat, 1)
#' mod2 <- mirt(dat, 1, itemtype = c(rep('2PL', 4), '3PL'))
#' anova(mod, mod2) #unbounded test
#' anova(mod, mod2, bounded = TRUE) #bounded
#'
#' # priors
#' model <- 'F = 1-5
#' PRIOR = (5, g, norm, -1, 1)'
#' mod1b <- mirt(dat, model, itemtype = c(rep('2PL', 4), '3PL'))
#' anova(mod1b)
#'
#' model2 <- 'F = 1-5
#' PRIOR = (1-5, g, norm, -1, 1)'
#' mod2b <- mirt(dat, model2, itemtype = '3PL')
#' anova(mod1b, mod2b)
#'
#' }
setMethod(
f = "anova",
signature = signature(object = 'SingleGroupClass'),
definition = function(object, object2, ...,
bounded = FALSE, mix = 0.5, frame = 1, verbose = FALSE){
if(frame > 1){
nms1 <- deparse(substitute(object, env = parent.frame(frame+1)))
nms2 <- deparse(substitute(object2, env = parent.frame(frame)))
} else {
nms1 <- deparse(substitute(object, env = parent.frame()))
nms2 <- deparse(substitute(object2, env = environment()))
}
dots <- list(...)
if(length(dots)){
nms3 <- if(frame > 1)
deparse(substitute(list(...), env = parent.frame(frame)))
else deparse(substitute(list(...), env = environment()))
nms3 <- gsub("list\\(", "", nms3)
nms3 <- gsub(")", "", nms3)
nms3 <- strsplit(nms3, ", ")[[1]]
dots <- c(object, object2, dots)
ret <- vector('list', length(dots)-1L)
for(i in 1L:length(ret)){
ret[[i]] <- anova(dots[[i]], dots[[i+1L]], bounded=bounded,
mix=mix)
if(i > 1L)
ret[[i]] <- ret[[i]][2L, ]
}
ret <- do.call(rbind, ret)
rownames(ret) <- c(nms1, nms2, nms3)
return(ret)
}
if(missing(object2)){
hasPriors <- object@Fit$logPrior != 0
ret <- data.frame(AIC = object@Fit$AIC,
SABIC = object@Fit$SABIC,
HQ = object@Fit$HQ,
BIC = object@Fit$BIC,
logLik = object@Fit$logLik)
if(hasPriors)
ret$logPost = object@Fit$logPrior + object@Fit$logLik
ret <- as.mirt_df(ret)
rownames(ret) <- nms1
return(ret)
}
df <- object@Fit$df - object2@Fit$df
ret <- data.frame(AIC = c(object@Fit$AIC, object2@Fit$AIC),
SABIC = c(object@Fit$SABIC, object2@Fit$SABIC),
HQ = c(object@Fit$HQ, object2@Fit$HQ),
BIC = c(object@Fit$BIC, object2@Fit$BIC),
logLik = c(object@Fit$logLik, object2@Fit$logLik))
if(any(object2@Fit$logPrior != 0 || object@Fit$logPrior != 0)){
ret$logPost = c(object@Fit$logLik + object@Fit$logPrior,
object2@Fit$logLik + object2@Fit$logPrior)
ret$df <- c(NA, abs(df))
} else {
X2 <- 2*object2@Fit$logLik - 2*object@Fit$logLik
ret$X2 <- c(NA, X2)
ret$df <- c(NA, df)
ret$p <- c(NA, 1 - pchisq(X2,abs(df)))
if(bounded)
ret$p[2L] <- 1 - mixX2(X2, df=abs(df), mix=mix)
ret$p[ret$X2 < 0] <- NaN
ret$p[ret$df <= 0] <- NaN
}
rownames(ret) <- c(nms1, nms2)
ret <- as.mirt_df(ret)
ret
}
)
#' Compute model residuals
#'
#' Return model implied residuals for linear dependencies between items or at the person level.
#' If the latent trait density was approximated (e.g., Davidian curves, Empirical histograms, etc)
#' then passing \code{use_dentype_estimate = TRUE} will use the internally saved quadrature and
#' density components (where applicable).
#'
#' @param object an object of class \code{SingleGroupClass} or
#' \code{MultipleGroupClass}. Bifactor models are automatically detected and utilized for
#' better accuracy
#' @param type type of residuals to be displayed.
#' Can be either \code{'LD'} or \code{'LDG2'} for a local dependence matrix based on the
#' X2 or G2 statistics (Chen & Thissen, 1997), \code{'Q3'} for the statistic proposed by
#' Yen (1984), \code{'JSI'} for the jack-knife statistic proposed Edwards et al. (2018),
#' \code{'exp'} for the expected values for the frequencies of every response pattern,
#' and \code{'expfull'} for the expected values for every theoretically observable response pattern.
#' For the 'LD' and 'LDG2' types, the upper diagonal elements represent the standardized
#' residuals in the form of signed Cramers V coefficients
#' @param tables logical; for LD type, return the observed, expected, and standardized residual
#' tables for each item combination?
#' @param df.p logical; print the degrees of freedom and p-values?
#' @param approx.z logical; transform \eqn{\chi^2(df)} information from LD tests into approximate
#' z-ratios instead using the transformation \eqn{z=\sqrt{2 * \chi^2} - \sqrt{2 * df - 1}}?
#' @param full.scores logical; compute relevant statistics
#' for each subject in the original data?
#' @param printvalue a numeric value to be specified when using the \code{res='exp'}
#' option. Only prints patterns that have standardized residuals greater than
#' \code{abs(printvalue)}. The default (NULL) prints all response patterns
#' @param verbose logical; allow information to be printed to the console?
#' @param Theta a matrix of factor scores used for statistics that require empirical estimates (i.e., Q3).
#' If supplied, arguments typically passed to \code{fscores()} will be ignored and these values will
#' be used instead
#' @param theta_lim range for the integration grid
#' @param p.adjust method to use for adjusting all p-values (see \code{\link{p.adjust}}
#' for available options). Default is \code{'none'}
#' @param fold logical; apply the sum 'folding' described by Edwards et al. (2018) for the JSI statistic?
#' @param quadpts number of quadrature nodes to use. The default is extracted from model (if available)
#' or generated automatically if not available
#' @param QMC logical; use quasi-Monte Carlo integration? If \code{quadpts} is omitted the
#' default number of nodes is 5000
#' @param suppress a numeric value indicating which parameter local dependency combinations
#' to flag as being too high (for LD, LDG2, and Q3 the standardize correlations are used; for
#' JSI, the z-ratios are used). Absolute values for the standardized estimates greater than
#' this value will be returned, while all values less than this value will be set to missing
#' @param upper logical; which portion of the matrix (upper versus lower triangle)
#' should the \code{suppress} argument be applied to?
#' @param technical list of technical arguments when models are re-estimated (see \code{\link{mirt}}
#' for details)
#' @param ... additional arguments to be passed to \code{fscores()}
#'
#' @name residuals-method
#' @aliases residuals,SingleGroupClass-method residuals,MixtureClass-method
#' residuals,MultipleGroupClass-method residuals,DiscreteClass-method
#' @rdname residuals-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#'
#' Chen, W. H. & Thissen, D. (1997). Local dependence indices for item pairs using item
#' response theory. \emph{Journal of Educational and Behavioral Statistics, 22}, 265-289.
#'
#' Edwards, M. C., Houts, C. R. & Cai, L. (2018). A Diagnostic Procedure to Detect Departures
#' From Local Independence in Item Response Theory Models.
#' \emph{Psychological Methods, 23}, 138-149.
#'
#' Yen, W. (1984). Effects of local item dependence on the fit and equating performance of the three
#' parameter logistic model. \emph{Applied Psychological Measurement, 8}, 125-145.
#' @examples
#'
#' \dontrun{
#'
#' x <- mirt(Science, 1)
#' residuals(x)
#' residuals(x, tables = TRUE)
#' residuals(x, type = 'exp')
#' residuals(x, suppress = .15)
#' residuals(x, df.p = TRUE)
#' residuals(x, df.p = TRUE, p.adjust = 'fdr') # apply FWE control
#'
#' # Pearson's X2 estimate for goodness-of-fit
#' full_table <- residuals(x, type = 'expfull')
#' head(full_table)
#' X2 <- with(full_table, sum((freq - exp)^2 / exp))
#' df <- nrow(full_table) - extract.mirt(x, 'nest') - 1
#' p <- pchisq(X2, df = df, lower.tail=FALSE)
#' data.frame(X2, df, p, row.names='Pearson-X2')
#'
#' # above FOG test as a function
#' PearsonX2 <- function(x){
#' full_table <- residuals(x, type = 'expfull')
#' X2 <- with(full_table, sum((freq - exp)^2 / exp))
#' df <- nrow(full_table) - extract.mirt(x, 'nest') - 1
#' p <- pchisq(X2, df = df, lower.tail=FALSE)
#' data.frame(X2, df, p, row.names='Pearson-X2')
#' }
#' PearsonX2(x)
#'
#'
#' # extract results manually
#' out <- residuals(x, df.p = TRUE, verbose=FALSE)
#' str(out)
#' out$df.p[1,2]
#'
#' # with and without supplied factor scores
#' Theta <- fscores(x)
#' residuals(x, type = 'Q3', Theta=Theta)
#' residuals(x, type = 'Q3', method = 'ML')
#'
#' # Edwards et al. (2018) JSI statistic
#' N <- 250
#' a <- rnorm(10, 1.7, 0.3)
#' d <- rnorm(10)
#' dat <- simdata(a, d, N=250, itemtype = '2PL')
#'
#' mod <- mirt(dat, 1)
#' residuals(mod, type = 'JSI')
#' residuals(mod, type = 'JSI', fold=FALSE) # unfolded
#'
#' # LD between items 1-2
#' aLD <- numeric(10)
#' aLD[1:2] <- rnorm(2, 2.55, 0.15)
#' a2 <- cbind(a, aLD)
#' dat <- simdata(a2, d, N=250, itemtype = '2PL')
#'
#' mod <- mirt(dat, 1)
#'
#' # JSI executed in parallel over multiple cores
#' if(interactive()) mirtCluster()
#' residuals(mod, type = 'JSI')
#'
#' }
setMethod(
f = "residuals",
signature = signature(object = 'SingleGroupClass'),
definition = function(object, type = 'LD', p.adjust = 'none',
df.p = FALSE, approx.z = FALSE,
full.scores = FALSE, QMC = FALSE,
printvalue = NULL, tables = FALSE, verbose = TRUE, Theta = NULL,
suppress = NA, theta_lim = c(-6, 6), quadpts = NULL, fold = TRUE,
upper = TRUE, technical = list(), ...)
{
dots <- list(...)
if(.hasSlot(object@Model$lrPars, 'beta'))
stop('Latent regression models not yet supported')
discrete <- use_dentype_estimate <- FALSE
if(!is.null(dots$use_dentype_estimate))
use_dentype_estimate <- dots$use_dentype_estimate
if(!is.null(dots$discrete) || use_dentype_estimate) discrete <- TRUE
K <- object@Data$K
data <- object@Data$data
N <- nrow(data)
J <- ncol(data)
nfact <- ncol(object@Fit$F)
res <- matrix(0,J,J)
diag(res) <- NA
colnames(res) <- rownames(res) <- colnames(data)
if(!is.null(Theta)){
if(!is.matrix(Theta)) Theta <- matrix(Theta)
if(nrow(Theta) > nrow(data))
Theta <- Theta[-extract.mirt(object, 'completely_missing'), , drop=FALSE]
stopifnot("Theta does not have the correct number of rows" =
nrow(Theta) == nrow(object@Data$data))
}
if(!discrete){
if(is.null(quadpts)){
if(QMC) quadpts <- 5000L
else quadpts <- object@Options$quadpts
}
if(is.nan(quadpts))
quadpts <- select_quadpts(nfact)
theta <- as.matrix(seq(theta_lim[1L], theta_lim[2L], length.out = quadpts))
if(type != 'Q3'){
Theta <- if(QMC)
QMC_quad(npts=quadpts, nfact=nfact, lim=theta_lim)
else {
if(nfact > 3L)
warning('High-dimensional models should use QMC integration instead',
call.=FALSE)
thetaComb(theta, nfact)
}
} else if(is.null(Theta)){
Theta <- fscores(object, verbose=FALSE, full.scores=TRUE, leave_missing=TRUE, ...)
}
} else {
Theta <- object@Model$Theta
}
itemnames <- colnames(data)
listtabs <- list()
calcG2 <- ifelse(type == 'LDG2', TRUE, FALSE)
if(type %in% c('LD', 'LDG2')){
if(!discrete){
groupPars <- ExtractGroupPars(object@ParObjects$pars[[object@Data$nitems + 1L]])
if(QMC){
Theta <- Theta_meanSigma_shift(Theta, groupPars$gmeans, groupPars$gcov)
prior <- rep(1/nrow(Theta), nrow(Theta))
} else {
prior <- mirt_dmvnorm(Theta,groupPars$gmeans, groupPars$gcov)
prior <- prior/sum(prior)
}
} else {
prior <- object@Internals$Prior[[1L]]
}
df <- (object@Data$K - 1) %o% (object@Data$K - 1)
diag(df) <- NA
colnames(df) <- rownames(df) <- colnames(res)
for(i in seq_len(J)){
for(j in seq_len(J)){
if(i < j){
P1 <- ProbTrace(x=object@ParObjects$pars[[i]], Theta=Theta)
P2 <- ProbTrace(x=object@ParObjects$pars[[j]], Theta=Theta)
tab <- table(data[,i], data[,j], useNA = 'no')
Etab <- matrix(0,K[i],K[j])
NN <- sum(tab)
for(k in seq_len(K[i]))
for(m in seq_len(K[j]))
Etab[k,m] <- NN * sum(P1[,k] * P2[,m] * prior)
s <- try(gamma_cor(tab) - gamma_cor(Etab), TRUE)
if(is.nan(s) || is(s, 'try-error')){
res[i,j] <- res[j,i] <- NaN
next
}
if(s == 0) s <- 1
if(calcG2){
tmp <- tab
tmp[tab == 0] <- NA
res[j,i] <- 2 * sum(tmp * log(tmp/Etab), na.rm=TRUE)
} else {
res[j,i] <- sum(((tab - Etab)^2)/Etab)
}
res[i,j] <- sign(s) * sqrt( abs(res[j,i]) / (NN * min(c(K[i],K[j]) - 1L)))
df[i,j] <- pchisq(abs(res[j,i]), df=df[j,i], lower.tail=FALSE)
if(tables){
tmp <- paste0(itemnames[i], '_', itemnames[j])
listtabs[[tmp]] <- list(Obs=tab, Exp=Etab, std_res=(tab-Etab)/sqrt(Etab))
}
if(approx.z)
res[j,i] <- sign(res[i,j]) * (sqrt(2 * res[j,i]) - sqrt(2 * df[j,i] - 1))
}
}
}
if(df.p)
df[upper.tri(df)] <- p.adjust(df[upper.tri(df)], method=p.adjust)
if(tables) return(listtabs)
if(df.p){
class(df) <- c('mirt_matrix', 'matrix')
if(verbose){
cat("Degrees of freedom (lower triangle) and p-values:\n\n")
print(df, ..., na.print = " ")
cat("\n")
}
}
if(verbose){
cat("LD matrix (lower triangle) and standardized values.\n")
cat("\nUpper triangle summary:\n")
print(round(summary(res[upper.tri(res)]), 3))
cat("\n")
}
res <- suppressMat(res, suppress=suppress, upper=upper)
class(res) <- c('mirt_matrix', 'matrix')
if(verbose) print(res, ..., na.print = " ")
if(df.p){
ret <- list(df, res)
names(ret) <- c('df.p', type)
return(invisible(ret))
}
return(invisible(res))
} else if(type == 'exp'){
r <- object@Data$Freq[[1L]]
res <- (r - object@Internals$Pl * nrow(object@Data$data)) /
sqrt(object@Internals$Pl * nrow(object@Data$data))
expected <- N * object@Internals$Pl
tabdata <- object@Data$tabdata
ISNA <- is.na(rowSums(tabdata))
expected[ISNA] <- res[ISNA] <- NA
tabdata <- data.frame(tabdata,object@Data$Freq[[1L]],expected,res)
colnames(tabdata) <- c(colnames(object@Data$tabdata),"freq","exp","std.res")
if(full.scores){
tabdata[, 'exp'] <- object@Internals$Pl / r * N
tabdata2 <- object@Data$tabdata
stabdata2 <- apply(tabdata2, 1, paste, sep='', collapse = '/')
sfulldata <- apply(object@Data$data, 1, paste, sep='', collapse = '/')
scoremat <- tabdata[match(sfulldata, stabdata2), 'exp', drop = FALSE]
res <- (1-scoremat) / sqrt(scoremat)
colnames(res) <- 'std.res'
ret <- cbind(object@Data$data, scoremat, res)
ret[is.na(rowSums(ret)), c('exp', 'std.res')] <- NA
rownames(ret) <- NULL
ret <- as.mirt_df(ret)
return(ret)
} else {
tabdata <- tabdata[do.call(order, as.data.frame(tabdata[,1:J])),]
if(!is.null(printvalue)){
if(!is.numeric(printvalue)) stop('printvalue is not a number.', call.=FALSE)
tabdata <- tabdata[abs(tabdata[ ,ncol(tabdata)]) > printvalue, ]
}
tabdata <- as.mirt_df(tabdata)
rownames(tabdata) <- NULL
return(tabdata)
}
} else if(type == 'expfull'){
K <- extract.mirt(object, 'K')
nitems <- length(K)
resp <- vector('list', nitems)
for(i in seq_len(nitems))
resp[[i]] <- 0L:(K[i]-1L)
tabdata <- expand.grid(resp)
rownames(tabdata) <- NULL
tabdata <- t(t(tabdata) + extract.mirt(object, 'mins'))
colnames(tabdata) <- colnames(extract.mirt(object, 'data'))
sv <- mod2values(object)
itemtype <- extract.mirt(object, 'itemtype')
nfact <- extract.mirt(object, 'nfact')
tmpdat <- matrix(0, nrow=2, ncol=nitems)
colnames(tmpdat) <- colnames(tabdata)
large <- mirt(tmpdat, nfact, itemtype=itemtype, pars=sv, TOL=NaN, large='return')
large$tabdata <- poly2dich(tabdata)
large$Freq$all <- rep(1L, nrow(tabdata))
large$tabdata2 <- matrix(1L)
full_object <- mirt(tabdata, nfact, itemtype=itemtype, pars=sv, TOL=NaN, large=large)
Pl <- full_object@Internals$Pl
r <- integer(length(Pl))
ro <- object@Data$Freq[[1L]]
N <- sum(ro)
for(i in seq_len(length(ro))){
pick <- colSums(t(tabdata) == object@Data$tabdata[i,]) == ncol(tabdata)
r[pick] <- ro[i]
}
res <- (r - Pl * N) / sqrt(Pl * N)
expected <- N * Pl
tabdata <- data.frame(tabdata,r,expected,res)
colnames(tabdata) <- c(colnames(object@Data$data),"freq","exp","res")
tabdata <- tabdata[do.call(order, as.data.frame(tabdata[,1:J])),]
rownames(tabdata) <- 1:length(r)
if(!is.null(printvalue)){
if(!is.numeric(printvalue)) stop('printvalue is not a number.', call.=FALSE)
tabdata <- tabdata[abs(tabdata[ ,ncol(tabdata)]) > printvalue, ]
}
tabdata <- as.mirt_df(tabdata)
return(tabdata)
} else if(type == 'Q3'){
if(discrete && !use_dentype_estimate)
stop('residual type not supported for discrete density forms', call.=FALSE)
dat <- matrix(NA, N, 2L)
diag(res) <- 1
for(i in seq_len(J)){
ei <- extract.item(object, item=i)
EI <- expected.item(ei, Theta=Theta)
dat[ ,1L] <- object@Data$data[ ,i] - EI
for(j in seq_len(J)){
if(i < j){
ej <- extract.item(object, item=j)
EJ <- expected.item(ej, Theta=Theta)
dat[,2L] <- object@Data$data[ ,j] - EJ
tmpdat <- na.omit(dat)
res[i,j] <- res[j,i] <- cor(tmpdat)[1L,2L]
}
}
}
if(verbose){
cat("Q3 summary statistics:\n")
print(round(summary(res[upper.tri(res)]), 3))
cat("\n")
}
res <- suppressMat(res, suppress=suppress, upper=upper)
class(res) <- c('mirt_matrix', 'matrix')
if(verbose) print(res, ..., na.print = " ")
return(invisible(res))
} else if(type == 'JSI'){
nfact <- extract.mirt(object, 'nfact')
stopifnot(nfact == 1L)
nitems <- extract.mirt(object, 'nitems')
technical$omp <- FALSE
as_drop <- myLapply(seq_len(nitems), function(item, mod, technical, ...){
itemtype <- extract.mirt(mod, 'itemtype')[-item]
tmpdat <- extract.mirt(mod, 'data')[,-item]
tmpmod <- mirt(tmpdat, 1L, itemtype=itemtype, SE=TRUE, verbose=FALSE,
technical=technical, ...)
ret <- sapply(coef(tmpmod, printSE=TRUE)[1:ncol(tmpdat)], function(x) x[1L:2L, 'a1'])
ret
}, mod=object, technical=technical, ...)
as <- sapply(coef(object)[1:nitems], function(x) x[1L, 'a1'])
retmat <- matrix(NA, nitems, nitems)
colnames(retmat) <- rownames(retmat) <- extract.mirt(object, 'itemnames')
for(i in seq_len(nitems)){
tmp <- as_drop[[i]]
pick <- colnames(tmp)
zs <- (as[pick] - tmp[1L, ]) / tmp[2L, ]
retmat[i, pick] <- zs
}
if(fold) retmat <- retmat + t(retmat)
retmat <- suppressMat(retmat, suppress=suppress, upper=upper)
class(retmat) <- c('mirt_matrix', 'matrix')
if(verbose){
cat("JSI summary statistics:\n")
print(round(summary(na.omit(as.vector(retmat))), 3))
cat("\n")
print(retmat, ..., na.print = " ")
}
return(invisible(retmat))
} else {
stop('specified type does not exist', call.=FALSE)
}
}
)
#' Plot various test-implied functions from models
#'
#' Plot various test implied response functions from models estimated in the mirt package.
#'
#' @param x an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{DiscreteClass}
#' @param y an arbitrary missing argument required for \code{R CMD check}
#' @param type type of plot to view. Can be
#' \describe{
#' \item{\code{'info'}}{test information function}
#' \item{\code{'rxx'}}{for the reliability function}
#' \item{\code{'infocontour'}}{for the test information contours}
#' \item{\code{'SE'}}{for the test standard error function}
#' \item{\code{'infotrace'}}{item information traceline plots}
#' \item{\code{'infoSE'}}{a combined test information and standard error plot}
#' \item{\code{'trace'}}{item probability traceline plots}
#' \item{\code{'itemscore'}}{item scoring traceline plots}
#' \item{\code{'score'}}{expected total score surface}
#' \item{\code{'scorecontour'}}{expected total score contour plot}
#' \item{\code{'posteriorTheta'}}{posterior for the latent trait distribution}
#' \item{\code{'EAPsum'}}{compares sum-scores to the expected values based
#' on the EAP for sum-scores method (see \code{\link{fscores}})}
#' }
#'
#' Note that if \code{dentype = 'empiricalhist'} was used in estimation then
#' the type \code{'empiricalhist'}
#' also will be available to generate the empirical histogram plot, and if
#' \code{dentype = 'Davidian-#'} was used then the type \code{'Davidian'}
#' will also be available to generate the curve estimates at the quadrature
#' nodes used during estimation
#' @param drop2 logical; where appropriate, for dichotomous response items drop the lowest category
#' and provide information pertaining only to the second response option?
#' @param degrees numeric value ranging from 0 to 90 used in \code{plot} to compute angle
#' for information-based plots with respect to the first dimension.
#' If a vector is used then a bubble plot is created with the summed information
#' across the angles specified (e.g., \code{degrees = seq(0, 90, by=10)})
#' @param theta_lim lower and upper limits of the latent trait (theta) to be evaluated, and is
#' used in conjunction with \code{npts}
#' @param npts number of quadrature points to be used for plotting features.
#' Larger values make plots look smoother
#' @param MI a single number indicating how many imputations to draw to form bootstrapped confidence
#' intervals for the selected test statistic. If greater than 0 a plot will be drawn with a shaded
#' region for the interval
#' @param CI a number from 0 to 1 indicating the confidence interval to select when MI input is
#' used. Default uses the 95\% confidence (CI = .95)
#' @param rot allows rotation of the 3D graphics
#' @param which.items numeric vector indicating which items to be used when plotting. Default is
#' to use all available items
#' @param facet_items logical; apply grid of plots across items? If \code{FALSE}, items will be
#' placed in one plot for each group
#' @param profile logical; provide a profile plot of response probabilities (objects returned from
#' \code{\link{mdirt}} only)
#' @param auto.key plotting argument passed to \code{\link[lattice]{lattice}}
#' @param par.strip.text plotting argument passed to \code{\link[lattice]{lattice}}
#' @param par.settings plotting argument passed to \code{\link[lattice]{lattice}}
#' @param ehist.cut a probability value indicating a threshold for excluding cases in empirical
#' histogram plots. Values larger than the default will include more points in the tails of the
#' plot, potentially squishing the 'meat' of the plot to take up less area than visually desired
#' @param main argument passed to lattice. Default generated automatically
#' @param drape logical argument passed to lattice. Default generated automatically
#' @param colorkey logical argument passed to lattice. Default generated automatically
#' @param add.ylab2 logical argument passed to lattice. Default generated automatically
#' @param ... additional arguments to be passed to lattice
#'
#' @name plot-method
#' @export
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @aliases plot,SingleGroupClass-method
#' plot,MultipleGroupClass-method plot,SingleGroupClass,missing-method
#' plot,DiscreteClass,missing-method plot,MixtureClass,missing-method
#' @rdname plot-method
#' @examples
#'
#' \dontrun{
#' x <- mirt(Science, 1, SE=TRUE)
#' plot(x)
#' plot(x, type = 'info')
#' plot(x, type = 'infotrace')
#' plot(x, type = 'infotrace', facet_items = FALSE)
#' plot(x, type = 'infoSE')
#' plot(x, type = 'rxx')
#' plot(x, type = 'posteriorTheta')
#'
#' # confidence interval plots when information matrix computed
#' plot(x)
#' plot(x, MI=100)
#' plot(x, type='info', MI=100)
#' plot(x, type='SE', MI=100)
#' plot(x, type='rxx', MI=100)
#'
#' # use the directlabels package to put labels on tracelines
#' library(directlabels)
#' plt <- plot(x, type = 'trace')
#' direct.label(plt, 'top.points')
#'
#' # additional modifications can be made via update().
#' # See ?update.trellis for further documentation
#' plt
#' update(plt, ylab = expression(Prob(theta)),
#' main = "Item Traceline Functions") # ylab/main changed
#'
#' set.seed(1234)
#' group <- sample(c('g1','g2'), nrow(Science), TRUE)
#' x2 <- multipleGroup(Science, 1, group)
#' plot(x2)
#' plot(x2, type = 'trace')
#' plot(x2, type = 'trace', which.items = 1:2)
#' plot(x2, type = 'itemscore', which.items = 1:2)
#' plot(x2, type = 'trace', which.items = 1, facet_items = FALSE) #facet by group
#' plot(x2, type = 'info')
#'
#' x3 <- mirt(Science, 2)
#' plot(x3, type = 'info')
#' plot(x3, type = 'SE', theta_lim = c(-3,3))
#'
#' }
setMethod(
f = "plot",
signature = signature(x = 'SingleGroupClass', y = 'missing'),
definition = function(x, y, type = 'score', npts = 200, drop2 = TRUE, degrees = 45,
theta_lim = c(-6,6), which.items = 1:extract.mirt(x, 'nitems'),
MI = 0, CI = .95, rot = list(xaxis = -70, yaxis = 30, zaxis = 10),
facet_items = TRUE, main = NULL,
drape = TRUE, colorkey = TRUE, ehist.cut = 1e-10, add.ylab2 = TRUE,
par.strip.text = list(cex = 0.7),
par.settings = list(strip.background = list(col = '#9ECAE1'),
strip.border = list(col = "black")),
auto.key = list(space = 'right', points=FALSE, lines=TRUE),
profile = FALSE, ...)
{
dots <- list(...)
if(!(type %in% c('info', 'SE', 'infoSE', 'rxx', 'trace', 'score', 'itemscore',
'infocontour', 'infotrace', 'scorecontour', 'empiricalhist', 'Davidian',
'EAPsum', 'posteriorTheta')))
stop('type supplied is not supported')
if (any(degrees > 90 | degrees < 0))
stop('Improper angle specified. Must be between 0 and 90.', call.=FALSE)
rot <- list(x = rot[[1]], y = rot[[2]], z = rot[[3]])
nfact <- x@Model$nfact
if(nfact > 3) stop("Can't plot high dimensional solutions.", call.=FALSE)
J <- x@Data$nitems
if(x@ParObjects$pars[[J + 1L]]@dentype == 'custom')
theta_lim <- x@Internals$theta_lim
theta <- seq(theta_lim[1L],theta_lim[2L],length.out=npts/(nfact^2))
ThetaFull <- Theta <- thetaComb(theta, nfact)
prodlist <- attr(x@ParObjects$pars, 'prodlist')
if(all(x@Data$K[which.items] == 2L) && facet_items) auto.key <- FALSE
if(length(prodlist))
ThetaFull <- prodterms(Theta,prodlist)
if(length(degrees) > ncol(ThetaFull)) type <- 'infoangle'
if(length(degrees) == 1L) degrees <- rep(degrees, ncol(ThetaFull))
info <- numeric(nrow(ThetaFull))
if(type %in% c('info', 'infocontour', 'rxx', 'SE', 'infoSE'))
info <- testinfo(x, ThetaFull, degrees = degrees, which.items=which.items)
if(type == 'infoangle'){
for(i in seq_len(length(degrees)))
info <- info + testinfo(x, ThetaFull, degrees = rep(degrees[i], ncol(ThetaFull)),
which.items=which.items)
}
mins <- x@Data$mins
maxs <- extract.mirt(x, 'K') + mins - 1
rotate <- if(is.null(dots$rotate)) 'none' else dots$rotate
if (x@Options$exploratory){
if(!is.null(dots$rotate)){
so <- summary(x, verbose=FALSE, digits=5, ...)
a <- rotateLambdas(so) * 1.702
for(i in 1:J)
x@ParObjects$pars[[i]]@par[1:nfact] <- a[i, ]
}
}
itemtrace <- computeItemtrace(x@ParObjects$pars, ThetaFull, x@Model$itemloc,
CUSTOM.IND=x@Internals$CUSTOM.IND)
score <- c()
for(i in 1:J)
score <- c(score, (0:(x@Data$K[i]-1) + mins[i]) * (i %in% which.items))
score <- matrix(score, nrow(itemtrace), ncol(itemtrace), byrow = TRUE)
plt <- data.frame(cbind(info,score=rowSums(score*itemtrace),Theta=Theta))
bundle <- length(which.items) != J
gp <- ExtractGroupPars(x@ParObjects$pars[[J+1]])
if(MI > 0L && nfact == 1L){
tmpx <- x
if(!x@Options$SE)
stop('Must compute an information matrix', call.=FALSE)
covB <- x@vcov
pre.ev <- eigen(covB)
names <- colnames(covB)
tmp <- lapply(names, function(x, split){
as.numeric(strsplit(x, split=split)[[1L]][-1L])
}, split='\\.')
imputenums <- do.call(c, tmp)
CIscore <- CIinfo <- CIrxx <- matrix(0, MI, length(plt$score))
for(i in seq_len(MI)){
while(TRUE){
tmp <- try(imputePars(pars=x@ParObjects$pars, pre.ev=pre.ev,
imputenums=imputenums, constrain=x@Model$constrain),
silent=TRUE)
if(!is(tmp, 'try-error')) break
}
tmpx@ParObjects$pars <- tmp
gp2 <- ExtractGroupPars(tmp[[J+1]])
itemtrace <- computeItemtrace(tmpx@ParObjects$pars, ThetaFull, x@Model$itemloc,
CUSTOM.IND=x@Internals$CUSTOM.IND)
tmpscore <- rowSums(score * itemtrace)
CIscore[i, ] <- tmpscore
CIinfo[i, ] <- testinfo(tmpx, ThetaFull)
CIrxx[i, ] <- CIinfo[i, ] / (CIinfo[i, ] + 1/gp2$gcov[1L,1L])
}
}
mins <- mins[which.items]
maxs <- maxs[which.items]
ybump <- (max(maxs) - min(mins))/15
ybump_full <- (sum(maxs) - sum(mins))/15
if(type == 'EAPsum'){
if(is.null(main))
main <- "Expected vs Observed Sum-Scores"
fs <- fscores(x, method = 'EAPsum', full.scores=FALSE, verbose=FALSE, ...)
plt <- with(fs, data.frame(Scores=Sum.Scores, y=c(observed, expected),
group = rep(c('observed', 'expected'), each=nrow(fs))))
return(xyplot(y~Scores, plt, type='l', main = main, group=plt$group,
auto.key=auto.key, xlab = expression(Sum-Score), ylab=expression(n),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
if(nfact == 3){
colnames(plt) <- c("info", "score", "Theta1", "Theta2", "Theta3")
plt$SE <- 1 / sqrt(plt$info)
if(type == 'infocontour'){
if(is.null(main)){
main <- paste("Test Information Contour")
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(contourplot(info ~ Theta1 * Theta2 | Theta3, data = plt,
main = main, xlab = expression(theta[1]),
ylab = expression(theta[2]),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'scorecontour'){
if(is.null(main)){
main <- paste("Expected Score Contour")
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(contourplot(score ~ Theta1 * Theta2 | Theta3, data = plt,
ylim=c(sum(mins)-ybump_full, sum(maxs)+ybump_full),
main = main, xlab = expression(theta[1]),
ylab = expression(theta[2]), ylim=c(sum(mins)-.1, sum(maxs)+.1),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'info'){
if(is.null(main)){
main <- "Test Information"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(wireframe(info ~ Theta1 + Theta2 | Theta3, data = plt, main = main,
zlab=expression(I(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'SEcontour'){
if(is.null(main)){
main <- "Test Standard Errors"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(contourplot(score ~ Theta1 * Theta2 | Theta3, data = plt,
main = main, xlab = expression(theta[1]),
ylab = expression(theta[2]),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'score'){
if(is.null(main)){
main <- if(bundle) "Expected Bundle Score" else "Expected Total Score"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(wireframe(score ~ Theta1 + Theta2 | Theta3, data = plt, main = main,
ylim=c(sum(mins)-ybump_full, sum(maxs)+ybump_full),
zlab=expression(Total(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings,
ylim=c(sum(mins)-.1, sum(maxs)+.1), ...))
} else if(type == 'SE'){
if(is.null(main)){
main <- "Test Standard Errors"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(wireframe(SE ~ Theta1 + Theta2 | Theta3, data = plt, main = main,
zlab=expression(SE(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
stop('plot type not supported for three dimensional model', call.=FALSE)
}
} else if(nfact == 2){
colnames(plt) <- c("info", "score", "Theta1", "Theta2")
plt$SE <- 1 / sqrt(plt$info)
if(type == 'infocontour'){
if(is.null(main)){
main <- paste("Test Information Contour")
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(contourplot(info ~ Theta1 * Theta2, data = plt,
main = main, xlab = expression(theta[1]),
ylab = expression(theta[2]),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'scorecontour'){
if(is.null(main)){
main <- paste("Expected Score Contour")
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(contourplot(score ~ Theta1 * Theta2, data = plt,
ylim=c(sum(mins)-ybump_full, sum(maxs)+ybump_full),
main = main, xlab = expression(theta[1]),
ylab = expression(theta[2]),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'info'){
if(is.null(main)){
main <- "Test Information"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(wireframe(info ~ Theta1 + Theta2, data = plt, main = main,
zlab=expression(I(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'trace'){
if(is.null(main))
main <- 'Item Probability Functions'
P <- vector('list', length(which.items))
names(P) <- colnames(x@Data$data)[which.items]
ind <- 1L
for(i in which.items){
tmp <- probtrace(extract.item(x, i), ThetaFull)
if(ncol(tmp) == 2L && drop2) tmp <- tmp[,2, drop=FALSE]
tmp2 <- data.frame(P=as.numeric(tmp), cat=gl(ncol(tmp), k=nrow(ThetaFull),
labels=paste0('P', seq_len(ncol(tmp)))))
P[[ind]] <- tmp2
ind <- ind + 1L
}
nrs <- sapply(P, nrow)
Pstack <- do.call(rbind, P)
names <- c()
for(i in seq_len(length(nrs)))
names <- c(names, rep(names(P)[i], nrs[i]))
plotobj <- data.frame(Pstack, item=names, Theta=ThetaFull)
plotobj$item <- factor(plotobj$item, levels = colnames(x@Data$data)[which.items])
return(wireframe(P ~ Theta.1 * Theta.2 | item, plotobj, zlim = c(-0.1,1.1), group=plotobj$cat,
zlab=expression(P(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'infotrace'){
if(is.null(main))
main <- 'Item Information'
I <- matrix(NA, nrow(Theta), J)
for(i in which.items)
I[,i] <- iteminfo(extract.item(x, i), ThetaFull, degrees=degrees)
I <- t(na.omit(t(I)))
items <- rep(colnames(x@Data$data)[which.items], each=nrow(Theta))
plotobj <- data.frame(I = as.numeric(I), Theta=ThetaFull, item=items)
plotobj$item <- factor(plotobj$item, levels = colnames(x@Data$data)[which.items])
return(wireframe(I ~ Theta.1 * Theta.2 | item, plotobj, group=plotobj$cat,
zlab=expression(P(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'itemscore'){
if(is.null(main))
main <- 'Expected Item Score'
S <- vector('list', length(which.items))
names(S) <- colnames(x@Data$data)[which.items]
ind <- 1L
for(i in which.items){
S[[ind]] <- expected.item(extract.item(x, i), ThetaFull, mins[i])
ind <- ind + 1L
}
Sstack <- do.call(c, S)
names <- rep(names(S), each = nrow(ThetaFull))
plotobj <- data.frame(S=Sstack, item=names, Theta=ThetaFull)
plotobj$item <- factor(plotobj$item, levels = colnames(x@Data$data)[which.items])
return(wireframe(S ~ Theta.1 * Theta.2 | item, plotobj, group=plotobj$cat,
zlab=expression(P(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'SEcontour'){
if(is.null(main)){
main <- "Test Standard Errors"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(contourplot(score ~ Theta1 * Theta2, data = plt,
main = main, xlab = expression(theta[1]),
ylab = expression(theta[2]),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'score'){
if(is.null(main)){
main <- if(bundle) "Expected Bundle Score" else "Expected Total Score"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(wireframe(score ~ Theta1 + Theta2, data = plt, main = main,
zlim=c(sum(mins)-ybump_full, sum(maxs)+ybump_full),
zlab=expression(Total(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'infoangle'){
if(is.null(main)){
main <- 'Information across different angles'
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
graphics::symbols(plt$Theta1, plt$Theta2, circles = sqrt(plt$info/pi), inches = .35, fg='white', bg='blue',
xlab = expression(theta[1]), ylab = expression(theta[2]),
main = main)
} else if(type == 'SE'){
if(is.null(main)){
main <- "Test Standard Errors"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(wireframe(SE ~ Theta1 + Theta2, data = plt, main = main,
zlab=expression(SE(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
stop('plot type not supported for two dimensional model', call.=FALSE)
}
} else {
colnames(plt) <- c("info", "score", "Theta")
plt$SE <- 1 / sqrt(plt$info)
plt$rxx <- plt$info / (plt$info + 1/gp$gcov[1L,1L])
if(MI > 0){
bs_range <- function(x, CI){
ss <- sort(x)
N <- length(ss)
ret <- c(upper = ss[ceiling(N * (1 - (1-CI)/2))],
middle = median(x),
lower = ss[floor(N * (1-CI)/2)])
ret
}
tmp <- apply(CIscore, 2, bs_range, CI=CI)
plt$CIscoreupper <- tmp['upper', ]
plt$CIscorelower <- tmp['lower', ]
tmp <- apply(CIinfo, 2, bs_range, CI=CI)
plt$CIinfoupper <- tmp['upper', ]
plt$CIinfolower <- tmp['lower', ]
plt$CISElower <- 1/sqrt(tmp['upper', ])
plt$CISEupper <- 1/sqrt(tmp['lower', ])
tmp <- apply(CIrxx, 2, bs_range, CI=CI)
plt$CIrxxupper <- tmp['upper', ]
plt$CIrxxlower <- tmp['lower', ]
}
if(type == 'info'){
if(is.null(main))
main <- 'Test Information'
if(MI > 0){
return(xyplot(info ~ Theta, data=plt,
upper=plt$CIinfoupper, lower=plt$CIinfolower,
panel = function(x, y, lower, upper, ...){
panel.polygon(c(x, rev(x)), c(upper, rev(lower)),
col="#E6E6E6", border = FALSE, ...)
panel.xyplot(x, y, type='l', lty=1,...)
},
main = main, ylim=c(min(plt$CIinfolower), max(plt$CIinfoupper)),
ylab = expression(I(theta)), xlab = expression(theta),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
return(xyplot(info~Theta, plt, type='l', main = main,
xlab = expression(theta), ylab=expression(I(theta)),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
} else if(type == 'rxx'){
if(is.null(main))
main <- 'Reliability'
if(MI > 0){
return(xyplot(rxx ~ Theta, data=plt,
upper=plt$CIrxxupper, lower=plt$CIrxxlower,
panel = function(x, y, lower, upper, ...){
panel.polygon(c(x, rev(x)), c(upper, rev(lower)),
col="#E6E6E6", border = FALSE, ...)
panel.xyplot(x, y, type='l', lty=1,...)
},
main = main, ylim=c(-0.1, 1.1),
ylab = expression(r[xx](theta)), xlab = expression(theta),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
return(xyplot(rxx~Theta, plt, type='l', main = main, ylim=c(-0.1, 1.1),
xlab = expression(theta), ylab=expression(r[xx](theta)),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
} else if(type == 'SE'){
if(is.null(main))
main <- 'Test Standard Errors'
if(MI > 0){
return(xyplot(SE ~ Theta, data=plt,
upper=plt$CISEupper, lower=plt$CISElower,
panel = function(x, y, lower, upper, ...){
panel.polygon(c(x, rev(x)), c(upper, rev(lower)),
col="#E6E6E6", border = FALSE, ...)
panel.xyplot(x, y, type='l', lty=1,...)
},
main = main, ylim=c(min(plt$CISElower), max(plt$CISEupper)),
ylab = expression(I(theta)), xlab = expression(theta),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
return(xyplot(SE~Theta, plt, type='l', main = main,
xlab = expression(theta), ylab=expression(SE(theta)),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
} else if(type == 'infoSE'){
if(is.null(main))
main <- 'Test Information and Standard Errors'
par.settings <- c(par.settings,
lattice::simpleTheme(col = c("#0080ff",'red'), lty = 1:2))
obj1 <- xyplot(info~Theta, plt, type='l', main = main,
xlab = expression(theta), ylab=expression(I(theta)),
par.strip.text=par.strip.text, par.settings=par.settings)
obj2 <- xyplot(SE~Theta, plt, type='l', ylab=expression(SE(theta)),
par.strip.text=par.strip.text, par.settings=par.settings)
if(requireNamespace("latticeExtra", quietly = TRUE)){
return(latticeExtra::doubleYScale(obj1, obj2, add.ylab2 = add.ylab2, ...))
} else {
stop('latticeExtra package is not available. Please install.', call.=FALSE)
}
} else if(type == 'trace'){
if(is.null(main))
main <- 'Item Probability Functions'
P <- vector('list', length(which.items))
names(P) <- colnames(x@Data$data)[which.items]
ind <- 1L
alltwocats <- all(extract.mirt(x, 'K')[which.items] == 2L)
for(i in which.items){
tmp <- probtrace(extract.item(x, i), ThetaFull)
if(ncol(tmp) == 2L && (facet_items || (!facet_items && alltwocats)) && drop2)
tmp <- tmp[,2, drop=FALSE]
tmp2 <- data.frame(P=as.numeric(tmp), cat=gl(ncol(tmp), k=nrow(Theta),
labels=paste0('P', seq_len(ncol(tmp)))))
P[[ind]] <- tmp2
ind <- ind + 1L
}
nrs <- sapply(P, nrow)
Pstack <- do.call(rbind, P)
names <- c()
for(i in seq_len(length(nrs)))
names <- c(names, rep(names(P)[i], nrs[i]))
plotobj <- data.frame(Pstack, item=names, Theta=Theta)
plotobj$item <- factor(plotobj$item, levels = colnames(x@Data$data)[which.items])
if(facet_items){
return(xyplot(P ~ Theta|item, plotobj, ylim = c(-0.1,1.1), groups = cat,
xlab = expression(theta), ylab = expression(P(theta)),
auto.key = auto.key, type = 'l', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
return(xyplot(P ~ Theta|cat, plotobj, groups=plotobj$item, ylim = c(-0.1,1.1),
xlab = expression(theta), ylab = expression(P(theta)),
auto.key = auto.key, type = 'l', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
} else if(type == 'itemscore'){
if(is.null(main))
main <- 'Expected Item Score'
S <- vector('list', length(which.items))
names(S) <- colnames(x@Data$data)[which.items]
ind <- 1L
for(i in which.items){
S[[ind]] <- expected.item(extract.item(x, i), ThetaFull, mins[i])
ind <- ind + 1L
}
Sstack <- do.call(c, S)
names <- rep(names(S), each = nrow(ThetaFull))
plotobj <- data.frame(S=Sstack, item=names, Theta=Theta)
plotobj$item <- factor(plotobj$item, levels = colnames(x@Data$data)[which.items])
if(facet_items){
return(xyplot(S ~ Theta|item, plotobj, ylim=c(min(mins)-ybump, max(maxs)+ybump),
xlab = expression(theta), ylab = expression(S(theta)),
auto.key = auto.key, type = 'l', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
return(xyplot(S ~ Theta, plotobj, groups=plotobj$item, ylim=c(min(mins)-.1, max(maxs)+.1),
xlab = expression(theta), ylab = expression(S(theta)),
auto.key = auto.key, type = 'l', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
} else if(type == 'infotrace'){
if(is.null(main))
main <- 'Item Information'
I <- matrix(NA, nrow(Theta), J)
for(i in which.items)
I[,i] <- iteminfo(extract.item(x, i), ThetaFull)
I <- t(na.omit(t(I)))
items <- rep(colnames(x@Data$data)[which.items], each=nrow(Theta))
plotobj <- data.frame(I = as.numeric(I), Theta=Theta, item=items)
plotobj$item <- factor(plotobj$item, levels = colnames(x@Data$data)[which.items])
if(facet_items){
return(xyplot(I ~ Theta|item, plotobj,
xlab = expression(theta), ylab = expression(I(theta)),
auto.key = auto.key, type = 'l', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
return(xyplot(I ~ Theta, plotobj, groups = plotobj$item,
xlab = expression(theta), ylab = expression(I(theta)),
auto.key = auto.key, type = 'l', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
} else if(type == 'score'){
if(is.null(main))
main <- if(bundle) "Expected Bundle Score" else "Expected Total Score"
if(MI > 0){
return(xyplot(score ~ Theta, data=plt,
ylim=c(sum(mins)-ybump_full, sum(maxs)+ybump_full),
upper=plt$CIscoreupper, lower=plt$CIscorelower,
panel = function(x, y, lower, upper, ...){
panel.polygon(c(x, rev(x)), c(upper, rev(lower)),
col="#E6E6E6", border = FALSE, ...)
panel.xyplot(x, y, type='l', lty=1,...)
},
main = main,
ylab = expression(T(theta)), xlab = expression(theta),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
return(xyplot(score ~ Theta, plt, ylim=c(sum(mins)-ybump_full, sum(maxs)+ybump_full),
xlab = expression(theta), ylab = expression(T(theta)),
type = 'l', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
} else if(type == 'empiricalhist'){
if(!(x@Options$dentype %in% c('EH', 'EHW')))
stop('Empirical histogram was not estimated for this object', call.=FALSE)
if(is.null(main))
main <- 'Empirical Histogram'
Prior <- x@Internals$Prior[[1L]]
Theta <- x@Model$Theta
# Prior <- Prior * nrow(x@Data$data)
cuts <- cut(Theta, floor(npts/2))
Prior <- do.call(c, lapply(split(Prior, cuts), mean))
Theta <- do.call(c, lapply(split(Theta, cuts), mean))
keep1 <- min(which(Prior > ehist.cut))
keep2 <- max(which(Prior > ehist.cut))
plt <- data.frame(Theta = Theta, Prior = Prior)
plt <- plt[keep1:keep2, , drop=FALSE]
return(xyplot(Prior ~ Theta, plt,
xlab = expression(theta), ylab = 'Density',
type = 'b', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'Davidian'){
if(x@Options$dentype != 'Davidian')
stop('Davidian curve was not estimated for this object', call.=FALSE)
if(is.null(main))
main <- 'Davidian Curve'
Prior <- x@Internals$Prior[[1L]]
Theta <- x@Model$Theta
# Prior <- Prior * nrow(x@Data$data)
plt <- data.frame(Theta = Theta, Prior = Prior)
return(xyplot(Prior ~ Theta, plt,
xlab = expression(theta), ylab = 'Density',
type = 'b', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'posteriorTheta'){
if(is.null(main))
main <- 'Posterior Distribution of Latent Trait'
plt <- extract.mirt(x, 'thetaPosterior')[[1]]
return(xyplot(posterior ~ Theta, plt,
xlab = expression(theta), ylab = 'Density',
type = 'b', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}else {
stop('plot not supported for unidimensional models', call.=FALSE)
}
}
}
)
mirt2traditional <- function(x, vcov, nfact){
cls <- class(x)
opar <- par <- x@par
which.a <- which(x@par[1L:nfact] != 0)
if(length(which.a) != 1L) return(x)
a.nms <- if(nfact == 1L) 'a' else paste0('a', 1L:nfact)
if(cls != 'GroupPars')
ncat <- x@ncat
if(cls == 'dich'){
fns <- vector('list', nfact + 3L)
fns[[nfact+1L]] <- function(par, index, opar){
if(index == (nfact + 1L)){
opar[c(which.a, nfact + 1L)] <- par
ret <- -opar[nfact + 1L]/opar[which.a]
}
ret
}
fns[[nfact+2L]] <- function(par, index, opar){
if(index == nfact + 2L)
ret <- plogis(par)
ret
}
fns[[nfact+3L]] <- function(par, index, opar){
if(index == nfact + 3L)
ret <- plogis(par)
ret
}
delta_index <- c(as.list(rep(NA, nfact)),
list(c(which.a, nfact + 1L),
nfact + 2L, nfact+3L))
par[nfact + 1L] <- -par[nfact + 1L]/par[which.a]
par[nfact + 2L] <- plogis(par[nfact + 2L])
par[nfact + 3L] <- plogis(par[nfact + 3L])
names(par) <- c(a.nms, 'b', 'g', 'u')
} else if(cls == 'fivePL'){
fns <- vector('list', nfact + 4L)
fns[[nfact+1L]] <- function(par, index, opar){
if(index == (nfact + 1L)){
opar[c(which.a, nfact + 1L)] <- par
ret <- -opar[nfact + 1L]/opar[which.a]
}
ret
}
fns[[nfact+2L]] <- function(par, index, opar){
if(index == nfact + 2L)
ret <- plogis(par)
ret
}
fns[[nfact+3L]] <- function(par, index, opar){
if(index == nfact + 3L)
ret <- plogis(par)
ret
}
fns[[nfact+4L]] <- function(par, index, opar){
if(index == nfact + 4L)
ret <- exp(par)
ret
}
delta_index <- c(as.list(rep(NA, nfact)),
list(c(which.a, nfact + 1L),
nfact + 2L, nfact+3L, nfact+4L))
par[nfact + 1L] <- -par[nfact + 1L]/par[which.a]
par[nfact + 2L] <- plogis(par[nfact + 2L])
par[nfact + 3L] <- plogis(par[nfact + 3L])
par[nfact + 4L] <- exp(par[nfact + 4L])
names(par) <- c(a.nms, 'b', 'g', 'u', 'S')
} else if(cls == 'graded'){
fns <- vector('list', ncat + nfact-1L)
for(i in 2L:ncat - 1L){
fns[[i + nfact]] <- function(par, index, opar){
if(index > nfact){
opar[c(which.a, index)] <- par
ret <- -opar[index]/opar[which.a]
}
ret
}
}
delta_index <- vector('list', ncat + nfact - 1L)
for(i in 1:nfact)
delta_index[[i]] <- NA
for(i in 2:ncat-1L){
par[i+nfact] <- -par[i+nfact]/par[which.a]
delta_index[[i+nfact]] <- c(which.a, i+nfact)
}
names(par) <- c(a.nms, paste0('b', 2:ncat-1L))
} else if(cls == 'gpcm'){
fns <- vector('list', ncat+nfact)
for(i in 2L:ncat-1L){
fns[[i+nfact]] <- function(par, index, opar){
if(index > nfact){
if(index == (nfact+1))
opar[c(which.a, ncat + nfact + 2L)] <- par
else opar[c(which.a, ncat + index + nfact - 1L,
ncat + index + nfact)] <- par
par <- opar
ds <- par[(nfact+1):length(par)]/par[which.a]
ds <- ds[-seq_len(ncat)]
newd <- numeric(length(ds)-1L)
for(i in 2:length(ds))
newd[i-1L] <- -(ds[i] - ds[i-1L])
ret <- c(par[1:nfact], newd)
ret <- ret[index]
}
ret
}
}
delta_index <- vector('list', ncat+nfact-1L)
for(i in 1:nfact)
delta_index[[i]] <- NA
ds <- par[(nfact+1):length(par)]/par[which.a]
ds <- ds[-seq_len(ncat)]
newd <- numeric(length(ds)-1L)
for(i in 2:length(ds))
newd[i-1L] <- -(ds[i] - ds[i-1L])
tmp <- rbind(1:(ncat-1), 2:ncat)
for(i in 2:length(ds)){
newd[i-1L] <- -(ds[i] - ds[i-1L])
delta_index[[i+nfact-1L]] <- c(which.a, tmp[,i-1L] + ncat + nfact)
}
delta_index[[nfact+1L]] <- c(which.a, ncat + nfact + 2L)
par <- c(x@par[1:nfact], newd)
names(par) <- c(a.nms, paste0('b', 1:length(newd)))
x@est <- x@est[c(1:nfact, (ncat+nfact+2L):length(x@est))]
} else if(cls == 'nominal'){
if(nfact > 1L) return(x)
fns <- vector('list', ncat*2)
for(i in 2L:length(par)-1L){
fns[[i]] <- function(par, index, opar){
if(index <= floor(length(opar)/2)){
as <- par[2:(ncat+1)] * par[1]
as <- as - mean(as)
ret <- as[index]
} else {
ds <- par
ds <- ds - mean(ds)
ret <- ds[index-ncat]
}
ret
}
}
delta_index <- vector('list', ncat*2)
for(i in 1L:(ncat)) delta_index[[i]] <- 1L:(ncat+1L)
for(i in (ncat+1L):(ncat*2)) delta_index[[i]] <- (ncat+2L):length(par)
as <- par[2:(ncat+1)] * par[1]
as <- as - mean(as)
ds <- par[(ncat+2):length(par)]
ds <- ds - mean(ds)
par <- c(as, ds)
names(par) <- c(paste0('a', 1:ncat), paste0('c', 1:ncat))
x@est <- rep(TRUE, ncat*2)
x@SEpar <- rep(as.numeric(NA), ncat*2)
} else if(cls == 'nestlogit'){
if(nfact > 1L) return(x)
fns <- vector('list', ncat*2 + 4)
fns[[2]] <- function(par, index, opar){
if(index == 2L){
opar[1L:2L] <- par
ret <- -opar[2L]/opar[1L]
}
ret
}
fns[[3]] <- function(par, index, opar){
if(index == 3L)
ret <- plogis(par)
ret
}
fns[[4]] <- function(par, index, opar){
if(index == 4L)
ret <- plogis(par)
ret
}
for(i in (2L:length(par)-1L) + 4){
fns[[i]] <- function(par, index, opar){
if(index <= (5 + ncat-2)){
as <- par
as <- as - mean(as)
ret <- as[index-4]
} else {
ds <- par
ds <- ds - mean(ds)
ret <- ds[index-4-(ncat-1)]
}
ret
}
}
delta_index <- vector('list', (ncat-1)*2)
for(i in 1L:(ncat-1)) delta_index[[i]] <- 1L:(ncat-1) + 4
for(i in ncat:((ncat-1)*2)) delta_index[[i]] <- 1L:(ncat-1) + 4 + ncat - 1L
delta_index <- c(list(NA, 1L:2L, 3L, 4L), delta_index)
par1 <- par[1:4]
par1[2] <- -par1[2]/par1[1]
par1[3] <- plogis(par1[3])
par1[4] <- plogis(par1[4])
names(par1) <- c(a.nms, 'b', 'g', 'u')
par2 <- par[5:length(par)]
as <- par2[1:(ncat-1)]
as <- as - mean(as)
ds <- par2[-c(1:(ncat-1))]
ds <- ds - mean(ds)
names(as) <- paste0('a', 1:(ncat-1))
names(ds) <- paste0('c', 1:(ncat-1))
par <- c(par1, as, ds)
x@est <- c(x@est[1:4], rep(TRUE, (ncat-1)*2))
x@SEpar <- c(x@SEpar[1], as.numeric(rep(NA, (ncat-1)*2 + 3)))
}
x@par <- par
names(x@est) <- names(par)
x@parnames <- names(x@par)
if(length(vcov) == 0L || (is.na(vcov[1L,1L]) || !(cls %in%
c('dich', 'graded', 'gpcm', 'nominal', 'nestlogit')))){
x@SEpar <- numeric()
} else {
nms <- colnames(vcov)
splt <- strsplit(nms, "\\.")
splt <- lapply(splt, function(x) as.integer(x[-1L]))
for(i in 1L:length(delta_index)){
if(!x@est[i]) next
if(is.na(delta_index[[i]][1L])) next
grad <- numerical_deriv(opar[delta_index[[i]]], fns[[i]], opar=opar, index=i)
parnum <- x@parnum[delta_index[[i]]]
pick <- numeric(length(grad))
for(j in 1L:length(parnum))
pick[j] <- suppressWarnings(min(which(do.call(c, lapply(splt, function(x)
any(x %in% parnum[j]))))))
grad <- grad[is.finite(pick)]
pick <- pick[is.finite(pick)]
x@SEpar[i] <- as.vector(sqrt(grad %*% vcov[pick, pick, drop=FALSE] %*% grad))
}
if(cls == 'gpcm') x@SEpar <- x@SEpar[1L:length(x@par)]
}
x
}
#' Convert traditional IRT metric into slope-intercept form used in mirt
#'
#' This is a helper function for users who have previously available traditional/classical
#' IRT parameters and want to know the equivalent slope-intercept translation used in \code{mirt}.
#' Note that this function assumes that the supplied models are unidimensional by definition (i.e.,
#' will have only one slope/discrimination) and in the logistic metric (i.e., logistic-ogive
#' scaling coefficient D=1). If there is no supported slope-intercept transformation
#' available then the original vector of parameters will be returned by default.
#'
#' Supported class transformations for the \code{cls} input are:
#'
#' \describe{
#' \item{Rasch, 2PL, 3PL, 3PLu, 4PL}{
#' Form must be: (discrimination, difficulty, lower-bound, upper-bound)
#' }
#' \item{graded}{
#' Form must be: (discrimination, difficulty 1, difficulty 2, ..., difficulty k-1)
#' }
#' \item{gpcm}{
#' Form must be: (discrimination, difficulty 1, difficulty 2, ..., difficulty k-1)
#' }
#' \item{nominal}{
#' Form must be: (discrimination 1, discrimination 2, ..., discrimination k,
#' difficulty 1, difficulty 2, ..., difficulty k)
#' }
#' }
#'
#' @param x a vector of parameters to transform
#' @param cls the class or itemtype of the supplied model
#' @param ncat the number of categories implied by the IRT model
#'
#' @return a named vector of slope-intercept parameters (if supported)
#' @export
#'
#' @examples
#'
#' # classical 3PL model
#' vec <- c(a=1.5, b=-1, g=.1, u=1)
#' slopeint <- traditional2mirt(vec, '3PL', ncat=2)
#' slopeint
#'
#' # classical graded model (four category)
#' vec <- c(a=1.5, b1=-1, b2=0, b3=1.5)
#' slopeint <- traditional2mirt(vec, 'graded', ncat=4)
#' slopeint
#'
#' # classical generalize partial credit model (four category)
#' vec <- c(a=1.5, b1=-1, b2=0, b3=1.5)
#' slopeint <- traditional2mirt(vec, 'gpcm', ncat=4)
#' slopeint
#'
#' # classical nominal model (4 category)
#' vec <- c(a1=.5, a2 = -1, a3=1, a4=-.5, d1=1, d2=-1, d3=-.5, d4=.5)
#' slopeint <- traditional2mirt(vec, 'nominal', ncat=4)
#' slopeint
#'
#'
traditional2mirt <- function(x, cls, ncat){
cls <- toInternalItemtype(cls)
if(cls == 'dich'){
par <- x
par[2L] <- -par[2L]*par[1L]
names(par) <- c('a1', 'd', 'g', 'u')
} else if(cls == 'graded'){
par <- x
for(i in 2L:ncat)
par[i] <- -par[i]*par[1L]
names(par) <- c('a1', paste0('d', 1:(length(par)-1)))
} else if(cls %in% c('gpcm', 'gpcmIRT')){
if(cls == 'gpcmIRT'){
x[-c(1, length(x))] <- x[-c(1, length(x))] - x[length(x)]
x <- x[-length(x)]
}
par <- c(x[1L], 0L:(ncat-1L), 0, x[-1L])
ds <- -par[-c(1:(ncat+1))]*par[1]
newd <- numeric(length(ds))
for(i in length(ds):2L)
newd[i] <- (ds[i] + ds[i-1L])
for(i in length(newd):3L){
pick <- if(i %% 2 == 0) seq(2, i, by = 2) else seq(1, i, by = 2)
newd[i] <- sum(newd[pick])
}
par <- c(par[1:(ncat+1)], newd)
names(par) <- c('a1', paste0('ak', 0:(ncat-1)), paste0('d', 0:(ncat-1)))
} else if(cls == 'nominal'){
as <- x[seq_len(length(x)/2)]
ds <- x[-seq_len(length(x)/2)]
a1 <- (as[ncat] - as[1L]) / (ncat-1L)
ak <- 1:ncat - 1
for(i in 2:(ncat-1))
ak[i] <- -(as[1L] - as[i]) / a1
dk <- ak
for(i in 2:ncat)
dk[i] <- ds[i] - ds[1L]
par <- c(a1, ak, dk)
names(par) <- c('a1', paste0('ak', 0:(ncat-1)), paste0('d', 0:(ncat-1)))
} else {
stop('traditional2mirt item class not supported', call.=FALSE)
}
par
}
#' Extract parameter variance covariance matrix
#'
#' Extract parameter variance covariance matrix
#'
#' @param object an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{MixedClass}
#'
#' @name vcov-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @export
#' @aliases vcov,SingleGroupClass-method vcov,MixtureClass-method
#' vcov,MultipleGroupClass-method vcov,MixedClass-method vcov,DiscreteClass-method
#' @rdname vcov-method
#' @examples
#'
#' \dontrun{
#' x <- mirt(Science, 1, SE=TRUE)
#' vcov(x)
#'
#' }
setMethod(
f = "vcov",
signature = signature(object = 'SingleGroupClass'),
definition = function(object){
extract.mirt(object, 'vcov')
}
)
#' Extract log-likelihood
#'
#' Extract the observed-data log-likelihood.
#'
#' @param object an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{MixedClass}
#'
#' @name logLik-method
#' @aliases logLik,SingleGroupClass-method logLik,MixtureClass-method
#' logLik,MultipleGroupClass-method logLik,MixedClass-method logLik,DiscreteClass-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @rdname logLik-method
#' @examples
#'
#' \dontrun{
#' x <- mirt(Science, 1)
#' logLik(x)
#'
#' }
setMethod(
f = "logLik",
signature = signature(object = 'SingleGroupClass'),
definition = function(object){
extract.mirt(object, 'logLik')
}
)
Try the mirt package in your browser
Any scripts or data that you put into this service are public.
mirt documentation built on Sept. 11, 2024, 7:14 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions traditional2mirt mirt2traditional
Documented in traditional2mirt
#' Print the model objects
#'
#' Print model object summaries to the console.
#'
#' @param x an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{MixedClass}
#'
#' @name print-method
#' @aliases print,SingleGroupClass-method
#' print,MultipleGroupClass-method print,MixedClass-method print,DiscreteClass-method
#' print,MixtureClass-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @rdname print-method
#' @examples
#'
#' \dontrun{
#' x <- mirt(Science, 1)
#' print(x)
#' }
setMethod(
f = "print",
signature = signature(x = 'SingleGroupClass'),
definition = function(x){
if(!length(x@time)){
cat('An object of class \"SingleGroupClass\"\n')
return(invisible(NULL))
}
cat("\nCall:\n", paste(deparse(x@Call), sep = "\n", collapse = "\n"),
"\n\n", sep = "")
cat("Full-information item factor analysis with ", x@Model$nfact, " factor(s).\n", sep="")
EMquad <- ''
if(x@Options$method == 'EM') EMquad <- c('\n using ', x@Options$quadpts, ' quadrature')
method <- x@Options$method
if(method == 'MIXED') method <- 'MHRM'
if(x@OptimInfo$converged)
cat("Converged within ", x@Options$TOL, ' tolerance after ', x@OptimInfo$iter, ' ',
method, " iterations.\n", sep = "")
else
cat("FAILED TO CONVERGE within ", x@Options$TOL, ' tolerance after ',
x@OptimInfo$iter, ' ', method, " iterations.\n", sep="")
cat('mirt version:', as.character(utils::packageVersion('mirt')), '\n')
cat('M-step optimizer:', x@Options$Moptim, '\n')
if(method %in% c('EM', 'QMCEM', 'BL', 'MCEM')){
if(method == 'EM' || method == 'QMCEM')
cat('EM acceleration:', x@Options$accelerate, '\n')
if(method == 'EM' || method == 'BL')
cat('Number of rectangular quadrature:', x@Options$quadpts)
else if(method == 'QMCEM')
cat('Number of quasi-Monte Carlo points:', x@Options$quadpts)
else if(method == 'MCEM')
cat('Number of Monte Carlo points:', x@Options$quadpts)
cat('\n')
}
dentype <- switch(x@Options$dentype,
"EH" = "Empirical histogram",
"EHW" = 'Empirical histogram (scaled)',
x@Options$dentype)
cat('Latent density type:', dentype, '\n')
if(method == 'MHRM')
cat("Average MH acceptance ratio(s):", paste0(round(x@OptimInfo$aveAR,3), collapse=', '), "\n")
if(!is.na(x@OptimInfo$secondordertest)){
cat("\nInformation matrix estimated with method:", x@Options$SE.type)
cat('\nSecond-order test: model ', if(!x@OptimInfo$secondordertest)
'is not a maximum or the information matrix is too inaccurate' else
'is a possible local maximum', sep = "")
if(x@OptimInfo$secondordertest)
cat("\nCondition number of information matrix = ", x@OptimInfo$condnum)
cat('\n')
}
if(length(x@Fit$logLik) > 0){
if(x@Fit$logPrior != 0){
cat("\nLog-posterior = ", x@Fit$logLik + x@Fit$logPrior, if(method == 'MHRM')
paste(', SE =', round(x@Fit$SElogLik,3)), "\n",sep='')
cat('Estimated parameters:', length(extract.mirt(x, 'parvec')), '\n')
} else {
cat("\nLog-likelihood = ", x@Fit$logLik, if(method == 'MHRM')
paste(', SE =', round(x@Fit$SElogLik,3)), "\n",sep='')
cat('Estimated parameters:', extract.mirt(x, 'nestpars'), '\n')
}
cat("AIC = ", x@Fit$AIC, "\n", sep='')
cat("BIC = ", x@Fit$BIC, "; SABIC = ", x@Fit$SABIC, "\n", sep='')
if(!is.nan(x@Fit$p)){
cat("G2 (", x@Fit$df,") = ", round(x@Fit$G2,2), ", p = ", round(x@Fit$p,4), sep='')
cat("\nRMSEA = ", round(x@Fit$RMSEA,3), ", CFI = ", round(x@Fit$CFI,3),
", TLI = ", round(x@Fit$TLI,3), sep='')
}
}
cat("\n")
}
)
#' Show model object
#'
#' Print model object summaries to the console.
#'
#' @param object an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{MixedClass}
#'
#' @name show-method
#' @aliases show,SingleGroupClass-method
#' show,MultipleGroupClass-method show,MixedClass-method show,DiscreteClass-method
#' show,MixtureClass-method
#' @rdname show-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @examples
#'
#' \dontrun{
#' x <- mirt(Science, 1)
#' show(x)
#' }
setMethod(
f = "show",
signature = signature(object = 'SingleGroupClass'),
definition = function(object){
print(object)
}
)
#' Summary of model object
#'
#' Transforms coefficients into a standardized factor loading's metric. For \code{MixedClass} objects,
#' the fixed and random coefficients are printed. Note that while the output to the console is rounded
#' to three digits, the returned list of objects is not. For simulations, use
#' \code{output <- summary(mod, verbose = FALSE)} to suppress the console messages.
#'
#' @param object an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{MixedClass}
#' @param rotate a string indicating which rotation to use for exploratory models, primarily
#' from the \code{GPArotation} package (see documentation therein).
#'
#' Rotations currently supported are: \code{'promax'}, \code{'oblimin'}, \code{'varimax'},
#' \code{'quartimin'}, \code{'targetT'}, \code{'targetQ'}, \code{'pstT'}, \code{'pstQ'},
#' \code{'oblimax'}, \code{'entropy'}, \code{'quartimax'}, \code{'simplimax'},
#' \code{'bentlerT'}, \code{'bentlerQ'}, \code{'tandemI'}, \code{'tandemII'},
#' \code{'geominT'}, \code{'geominQ'}, \code{'cfT'}, \code{'cfQ'}, \code{'infomaxT'},
#' \code{'infomaxQ'}, \code{'mccammon'}, \code{'bifactorT'}, \code{'bifactorQ'}.
#'
#' For models that are not exploratory this input will automatically be set to \code{'none'}
#' @param Target a dummy variable matrix indicting a target rotation pattern. This is required for
#' rotations such as \code{'targetT'}, \code{'targetQ'}, \code{'pstT'}, and \code{'pstQ'}
#' @param suppress a numeric value indicating which (possibly rotated) factor
#' loadings should be suppressed. Typical values are around .3 in most
#' statistical software. Default is 0 for no suppression
#' @param suppress.cor same as \code{suppress}, but for the correlation matrix
#' output
#' @param verbose logical; allow information to be printed to the console?
#' @param ... additional arguments to be passed
#'
#' @name summary-method
#' @aliases summary,SingleGroupClass-method
#' summary,MultipleGroupClass-method summary,MixedClass-method summary,DiscreteClass-method
#' summary,MixtureClass-method
#' @export
#' @rdname summary-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @seealso \code{\link{coef-method}}
#' @examples
#'
#' \dontrun{
#' x <- mirt(Science, 2)
#' summary(x)
#' summary(x, rotate = 'varimax')
#'
#' }
setMethod(
f = "summary",
signature = 'SingleGroupClass',
definition = function(object, rotate = 'oblimin', Target = NULL,
suppress = 0, suppress.cor = 0,
verbose = TRUE, ...){
if (!object@Options$exploratory || rotate == 'none') {
F <- object@Fit$F
F[abs(F) < suppress] <- NA
h2 <- as.matrix(object@Fit$h2)
SS <- apply(F^2,2,sum)
gp <- ExtractGroupPars(object@ParObjects$pars[[object@Data$nitems + 1L]])
Phi <- cov2cor(gp$gcov)
colnames(h2) <- "h2"
rownames(Phi) <- colnames(Phi) <- names(SS) <-
colnames(F)[seq_len(object@Model$nfact)]
loads <- cbind(F,h2)
if(verbose){
if(object@Options$exploratory)
cat("\nUnrotated factor loadings: \n\n")
print(loads, 3, na.print = " ")
cat("\nSS loadings: ", round(SS, 3), "\n")
cat("Proportion Var: ",round(SS/nrow(F), 3), "\n")
cat("\nFactor correlations: \n\n")
Phiprint <- Phi
Phiprint[abs(Phi) < suppress.cor] <- NA
Phiprint[upper.tri(Phiprint, diag = FALSE)] <- NA
print(round(Phiprint, 3), na.print = " ")
}
invisible(list(rotF=F,h2=h2,fcor=Phi))
} else {
F <- object@Fit$F
h2 <- as.matrix(object@Fit$h2)
colnames(h2) <- "h2"
rotF <- Rotate(F, rotate, Target = Target, ...)
SS <- apply(rotF$loadings^2,2,sum)
L <- rotF$loadings
L[abs(L) < suppress] <- NA
loads <- cbind(L,h2)
Phi <- diag(ncol(F))
if(!rotF$orthogonal){
Phi <- rotF$Phi
}
colnames(Phi) <- rownames(Phi) <- colnames(F)
if(verbose){
cat("\nRotation: ", rotate, "\n")
cat("\nRotated factor loadings: \n\n")
print(loads, 3, na.print = " ")
cat("\nRotated SS loadings: ",round(SS,3), "\n")
cat("\nFactor correlations: \n\n")
Phiprint <- Phi
Phiprint[abs(Phi) < suppress.cor] <- NA
Phiprint[upper.tri(Phiprint, diag = FALSE)] <- NA
print(round(Phiprint, 3), na.print = " ")
}
if(any(h2 > 1))
warning("Solution has Heywood cases. Interpret with caution.",
call.=FALSE)
invisible(list(rotF=rotF$loadings,h2=h2,fcor=Phi))
}
}
)
#' Extract raw coefs from model object
#'
#' Return a list (or data.frame) of raw item and group level coefficients. Note that while
#' the output to the console is rounded to three digits, the returned list of objects is not.
#' Hence, elements from \code{cfs <- coef(mod); cfs[[1]]} will contain the non-rounded results (useful
#' for simulations).
#'
#' @param object an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{MixedClass}
#' @param CI the amount of converged used to compute confidence intervals; default is
#' 95 percent confidence intervals
#' @param IRTpars logical; convert slope intercept parameters into traditional IRT parameters?
#' Only applicable to unidimensional models or models with simple structure (i.e., only one non-zero slope).
#' If a suitable ACOV estimate was computed in the fitted
#' model, and \code{printSE = FALSE}, then suitable CIs will be included based on the delta
#' method (where applicable)
#' @param rotate see \code{summary} method for details. The default rotation is \code{'none'}
#' @param Target a dummy variable matrix indicting a target rotation pattern
#' @param printSE logical; print the standard errors instead of the confidence intervals? When
#' \code{IRTpars = TRUE} then the delta method will be used to compute the associated standard errors
#' from mirt's default slope-intercept form
#' @param as.data.frame logical; convert list output to a data.frame instead?
#' @param simplify logical; if all items have the same parameter names (indicating they are
#' of the same class) then they are collapsed to a matrix, and a list of length 2 is returned
#' containing a matrix of item parameters and group-level estimates
#' @param unique return the vector of uniquely estimated parameters
#' @param verbose logical; allow information to be printed to the console?
#' @param rawug logical; return the untransformed internal g and u parameters?
#' If \code{FALSE}, g and u's are converted with the original format along with delta standard errors
#' @param ... additional arguments to be passed
#'
#' @name coef-method
#' @aliases coef,SingleGroupClass-method
#' coef,MultipleGroupClass-method coef,MixedClass-method coef,DiscreteClass-method
#' coef,MixtureClass-method
#' @rdname coef-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @export
#' @seealso \code{\link{summary-method}}
#' @examples
#'
#' \dontrun{
#' dat <- expand.table(LSAT7)
#' x <- mirt(dat, 1)
#' coef(x)
#' coef(x, IRTpars = TRUE)
#' coef(x, simplify = TRUE)
#'
#' #with computed information matrix
#' x <- mirt(dat, 1, SE = TRUE)
#' coef(x)
#' coef(x, printSE = TRUE)
#' coef(x, as.data.frame = TRUE)
#'
#' #two factors
#' x2 <- mirt(Science, 2)
#' coef(x2)
#' coef(x2, rotate = 'varimax')
#'
#' }
setMethod(
f = "coef",
signature = 'SingleGroupClass',
definition = function(object, CI = .95, printSE = FALSE, rotate = 'none', Target = NULL,
IRTpars = FALSE, rawug = FALSE, as.data.frame = FALSE,
simplify = FALSE, unique = FALSE, verbose = TRUE, ...){
dots <- list(...)
discrete <- ifelse(is.null(dots$discrete), FALSE, TRUE)
if(unique) return(extract.mirt(object, 'parvec'))
if(printSE && length(object@ParObjects$pars[[1L]]@SEpar)) rawug <- TRUE
if(CI >= 1 || CI <= 0)
stop('CI must be between 0 and 1', call.=FALSE)
z <- abs(qnorm((1 - CI)/2))
SEnames <- paste0('CI_', c((1 - CI)/2*100, ((1 - CI)/2 + CI)*100))
J <- object@Data$nitems
nfact <- object@Model$nfact + length(object@Model$prodlist)
a <- matrix(0, J, nfact)
for(i in 1:J)
a[i, ] <- ExtractLambdas(object@ParObjects$pars[[i]])
if (object@Options$exploratory && rotate != 'none'){
if(verbose) cat("\nRotation: ", rotate, "\n\n")
so <- summary(object, rotate=rotate, Target=Target, verbose=FALSE, ...)
a <- rotateLambdas(so) * 1.702
for(i in 1:J){
object@ParObjects$pars[[i]]@par[1:nfact] <- a[i, ]
object@ParObjects$pars[[i]]@SEpar <- numeric(0L)
}
object@ParObjects$pars[[J + 1]]@par[-c(1:nfact)] <- so$fcor[lower.tri(so$fcor, TRUE)]
}
if(IRTpars){
if(object@Model$nfact > 1L){
apars <- lapply(object@ParObjects$pars[-(extract.mirt(object, 'nitems') + 1L)],
function(x) x@par[1L:object@Model$nfact])
is_ss <- sapply(apars, function(x) sum(x != 0) == 1L)
if(!any(is_ss)){
warning(c('Traditional parameterization only available for unidimensional ',
'models or models with simple structure patterns (neither found)'), call.=FALSE)
IRTpars <- FALSE
}
}
if(IRTpars){
vcov <- vcov(object)
for(i in 1L:J){
if(class(object@ParObjects$pars[[i]]) %in% c('gpcmIRT')) next
object@ParObjects$pars[[i]] <- mirt2traditional(object@ParObjects$pars[[i]],
vcov=vcov, nfact=object@Model$nfact)
}
}
}
allPars <- list()
if(length(object@ParObjects$pars[[1L]]@SEpar) && !simplify){
if(printSE){
for(i in seq_len(J+1L)){
allPars[[i]] <- matrix(c(object@ParObjects$pars[[i]]@par,
object@ParObjects$pars[[i]]@SEpar),
2, byrow = TRUE)
rownames(allPars[[i]]) <- c('par', 'SE')
nms <- names(object@ParObjects$pars[[i]]@est)
if(i <= J && object@Model$itemtype[i] != 'custom' && !IRTpars){
nms[nms == 'g'] <- 'logit(g)'
nms[nms == 'u'] <- 'logit(u)'
}
colnames(allPars[[i]]) <- nms
}
} else {
for(i in seq_len(J+1L)){
allPars[[i]] <- matrix(c(object@ParObjects$pars[[i]]@par,
object@ParObjects$pars[[i]]@par -
z*object@ParObjects$pars[[i]]@SEpar,
object@ParObjects$pars[[i]]@par +
z*object@ParObjects$pars[[i]]@SEpar),
3, byrow = TRUE)
rownames(allPars[[i]]) <- c('par', SEnames)
colnames(allPars[[i]]) <- object@ParObjects$pars[[i]]@parnames
}
}
} else {
for(i in seq_len(J+1L)){
allPars[[i]] <- matrix(object@ParObjects$pars[[i]]@par, 1L)
colnames(allPars[[i]]) <- object@ParObjects$pars[[i]]@parnames
rownames(allPars[[i]]) <- 'par'
}
}
if(!rawug && !IRTpars){
allPars <- lapply(allPars, function(x, digits){
x[ , colnames(x) %in% c('g', 'u')] <-
antilogit(x[ , colnames(x) %in% c('g', 'u')])
x
})
}
names(allPars) <- c(colnames(object@Data$data), 'GroupPars')
if(as.data.frame)
allPars <- t(as.data.frame(allPars))
if(simplify && !as.data.frame){
allPars <- lapply(allPars, function(x) x[1L, , drop=FALSE])
items.old <- allPars[seq_len(length(allPars)-1L)]
nms <- lapply(items.old, colnames)
unms <- unique(do.call(c, nms))
items <- matrix(NA, length(items.old), length(unms))
rownames(items) <- names(items.old)
colnames(items) <- unms
for(i in seq_len(nrow(items)))
items[i, nms[[i]]] <- items.old[[i]]
nfact <- object@Model$nfact
if(discrete){
allPars <- list(items=items, group.intercepts=allPars$GroupPars)
} else if(object@ParObjects$pars[[J+1L]]@dentype != "custom"){
means <- allPars$GroupPars[seq_len(nfact)]
if(object@ParObjects$pars[[J+1L]]@dentype == "Davidian"){
covs <- matrix(NA, nfact, nfact)
covs[lower.tri(covs, TRUE)] <- allPars$GroupPars[2L]
covs[upper.tri(covs, FALSE)] <- covs[lower.tri(covs, FALSE)]
colnames(covs) <- rownames(covs) <- names(means) <-
object@Model$factorNames[seq_len(nfact)]
allPars <- list(items=items, means=means, cov=covs,
Davidian_phis=allPars$GroupPars[-c(1:2)])
} else {
covs <- matrix(NA, nfact, nfact)
if(object@ParObjects$pars[[J+1L]]@dentype == "mixture")
covs[lower.tri(covs, TRUE)] <-
allPars$GroupPars[-c(seq_len(nfact), length(allPars$GroupPars))]
else covs[lower.tri(covs, TRUE)] <- allPars$GroupPars[-seq_len(nfact)]
covs <- makeSymMat(covs)
colnames(covs) <- rownames(covs) <- names(means) <-
object@Model$factorNames[seq_len(nfact)]
allPars <- list(items=items, means=means, cov=covs)
}
} else allPars <- list(items=items, GroupPars=allPars$GroupPars)
}
if(.hasSlot(object@Model$lrPars, 'beta')){
allPars$lr.betas <- object@Model$lrPars@beta
if(!all(is.nan(object@Model$lrPars@SEpar))){
tmp <- allPars$lr.betas
if(printSE){
tmp[] <- object@Model$lrPars@SEpar
allPars$lr.betas <- list(betas=allPars$lr.betas, SE=tmp)
} else {
low <- tmp - z*object@Model$lrPars@SEpar
high <- tmp + z*object@Model$lrPars@SEpar
allPars$lr.betas <- list(betas=allPars$lr.betas, low, high)
names(allPars$lr.betas) <- c('betas', SEnames)
}
}
}
if(!as.data.frame)
class(allPars) <- c('mirt_list', 'list')
return(allPars)
}
)
#' Compare nested models with likelihood-based statistics
#'
#' Compare nested models using likelihood ratio test (X2), Akaike Information Criterion (AIC),
#' Bayesian Information Criterion (BIC),
#' Sample-Size Adjusted BIC (SABIC), and Hannan-Quinn (HQ) Criterion.
#' When given a sequence of objects, \code{anova} tests the models against one another
#' in the order specified. Note that the \code{object} inputs should be ordered in terms
#' of most constrained model to least constrained.
#'
#' @param object an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{MixedClass}, reflecting the most constrained model fitted
#' @param object2 a second model estimated from any of the mirt package estimation methods
#' @param ... additional less constrained model objects to be compared
#' sequentially to the previous model
#' @param bounded logical; are the two models comparing a bounded parameter (e.g., comparing a single
#' 2PL and 3PL model with 1 df)? If \code{TRUE} then a 50:50 mix of chi-squared distributions
#' is used to obtain the p-value
#' @param mix proportion of chi-squared mixtures. Default is 0.5
#' @param frame (internal parameter not for standard use)
#' @param verbose (deprecated argument)
#'
#' @return a \code{data.frame}/\code{mirt_df} object
#'
#' @name anova-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @aliases anova,SingleGroupClass-method
#' anova,MultipleGroupClass-method anova,MixedClass-method anova,DiscreteClass-method
#' anova,MixtureClass-method
#' @export
#' @rdname anova-method
#' @examples
#'
#' \dontrun{
#' x <- mirt(Science, 1)
#' x2 <- mirt(Science, 2)
#' anova(x, x2)
#'
#' # compare three models sequentially (X2 not always meaningful)
#' x3 <- mirt(Science, 1, 'gpcm')
#' x4 <- mirt(Science, 1, 'nominal')
#' anova(x, x2, x3, x4)
#'
#' # in isolation
#' anova(x)
#'
#' # with priors on first model
#' model <- "Theta = 1-4
#' PRIOR = (1-4, a1, lnorm, 0, 10)"
#' xp <- mirt(Science, model)
#' anova(xp, x2)
#' anova(xp)
#'
#' # bounded parameter
#' dat <- expand.table(LSAT7)
#' mod <- mirt(dat, 1)
#' mod2 <- mirt(dat, 1, itemtype = c(rep('2PL', 4), '3PL'))
#' anova(mod, mod2) #unbounded test
#' anova(mod, mod2, bounded = TRUE) #bounded
#'
#' # priors
#' model <- 'F = 1-5
#' PRIOR = (5, g, norm, -1, 1)'
#' mod1b <- mirt(dat, model, itemtype = c(rep('2PL', 4), '3PL'))
#' anova(mod1b)
#'
#' model2 <- 'F = 1-5
#' PRIOR = (1-5, g, norm, -1, 1)'
#' mod2b <- mirt(dat, model2, itemtype = '3PL')
#' anova(mod1b, mod2b)
#'
#' }
setMethod(
f = "anova",
signature = signature(object = 'SingleGroupClass'),
definition = function(object, object2, ...,
bounded = FALSE, mix = 0.5, frame = 1, verbose = FALSE){
if(frame > 1){
nms1 <- deparse(substitute(object, env = parent.frame(frame+1)))
nms2 <- deparse(substitute(object2, env = parent.frame(frame)))
} else {
nms1 <- deparse(substitute(object, env = parent.frame()))
nms2 <- deparse(substitute(object2, env = environment()))
}
dots <- list(...)
if(length(dots)){
nms3 <- if(frame > 1)
deparse(substitute(list(...), env = parent.frame(frame)))
else deparse(substitute(list(...), env = environment()))
nms3 <- gsub("list\\(", "", nms3)
nms3 <- gsub(")", "", nms3)
nms3 <- strsplit(nms3, ", ")[[1]]
dots <- c(object, object2, dots)
ret <- vector('list', length(dots)-1L)
for(i in 1L:length(ret)){
ret[[i]] <- anova(dots[[i]], dots[[i+1L]], bounded=bounded,
mix=mix)
if(i > 1L)
ret[[i]] <- ret[[i]][2L, ]
}
ret <- do.call(rbind, ret)
rownames(ret) <- c(nms1, nms2, nms3)
return(ret)
}
if(missing(object2)){
hasPriors <- object@Fit$logPrior != 0
ret <- data.frame(AIC = object@Fit$AIC,
SABIC = object@Fit$SABIC,
HQ = object@Fit$HQ,
BIC = object@Fit$BIC,
logLik = object@Fit$logLik)
if(hasPriors)
ret$logPost = object@Fit$logPrior + object@Fit$logLik
ret <- as.mirt_df(ret)
rownames(ret) <- nms1
return(ret)
}
df <- object@Fit$df - object2@Fit$df
ret <- data.frame(AIC = c(object@Fit$AIC, object2@Fit$AIC),
SABIC = c(object@Fit$SABIC, object2@Fit$SABIC),
HQ = c(object@Fit$HQ, object2@Fit$HQ),
BIC = c(object@Fit$BIC, object2@Fit$BIC),
logLik = c(object@Fit$logLik, object2@Fit$logLik))
if(any(object2@Fit$logPrior != 0 || object@Fit$logPrior != 0)){
ret$logPost = c(object@Fit$logLik + object@Fit$logPrior,
object2@Fit$logLik + object2@Fit$logPrior)
ret$df <- c(NA, abs(df))
} else {
X2 <- 2*object2@Fit$logLik - 2*object@Fit$logLik
ret$X2 <- c(NA, X2)
ret$df <- c(NA, df)
ret$p <- c(NA, 1 - pchisq(X2,abs(df)))
if(bounded)
ret$p[2L] <- 1 - mixX2(X2, df=abs(df), mix=mix)
ret$p[ret$X2 < 0] <- NaN
ret$p[ret$df <= 0] <- NaN
}
rownames(ret) <- c(nms1, nms2)
ret <- as.mirt_df(ret)
ret
}
)
#' Compute model residuals
#'
#' Return model implied residuals for linear dependencies between items or at the person level.
#' If the latent trait density was approximated (e.g., Davidian curves, Empirical histograms, etc)
#' then passing \code{use_dentype_estimate = TRUE} will use the internally saved quadrature and
#' density components (where applicable).
#'
#' @param object an object of class \code{SingleGroupClass} or
#' \code{MultipleGroupClass}. Bifactor models are automatically detected and utilized for
#' better accuracy
#' @param type type of residuals to be displayed.
#' Can be either \code{'LD'} or \code{'LDG2'} for a local dependence matrix based on the
#' X2 or G2 statistics (Chen & Thissen, 1997), \code{'Q3'} for the statistic proposed by
#' Yen (1984), \code{'JSI'} for the jack-knife statistic proposed Edwards et al. (2018),
#' \code{'exp'} for the expected values for the frequencies of every response pattern,
#' and \code{'expfull'} for the expected values for every theoretically observable response pattern.
#' For the 'LD' and 'LDG2' types, the upper diagonal elements represent the standardized
#' residuals in the form of signed Cramers V coefficients
#' @param tables logical; for LD type, return the observed, expected, and standardized residual
#' tables for each item combination?
#' @param df.p logical; print the degrees of freedom and p-values?
#' @param approx.z logical; transform \eqn{\chi^2(df)} information from LD tests into approximate
#' z-ratios instead using the transformation \eqn{z=\sqrt{2 * \chi^2} - \sqrt{2 * df - 1}}?
#' @param full.scores logical; compute relevant statistics
#' for each subject in the original data?
#' @param printvalue a numeric value to be specified when using the \code{res='exp'}
#' option. Only prints patterns that have standardized residuals greater than
#' \code{abs(printvalue)}. The default (NULL) prints all response patterns
#' @param verbose logical; allow information to be printed to the console?
#' @param Theta a matrix of factor scores used for statistics that require empirical estimates (i.e., Q3).
#' If supplied, arguments typically passed to \code{fscores()} will be ignored and these values will
#' be used instead
#' @param theta_lim range for the integration grid
#' @param p.adjust method to use for adjusting all p-values (see \code{\link{p.adjust}}
#' for available options). Default is \code{'none'}
#' @param fold logical; apply the sum 'folding' described by Edwards et al. (2018) for the JSI statistic?
#' @param quadpts number of quadrature nodes to use. The default is extracted from model (if available)
#' or generated automatically if not available
#' @param QMC logical; use quasi-Monte Carlo integration? If \code{quadpts} is omitted the
#' default number of nodes is 5000
#' @param suppress a numeric value indicating which parameter local dependency combinations
#' to flag as being too high (for LD, LDG2, and Q3 the standardize correlations are used; for
#' JSI, the z-ratios are used). Absolute values for the standardized estimates greater than
#' this value will be returned, while all values less than this value will be set to missing
#' @param upper logical; which portion of the matrix (upper versus lower triangle)
#' should the \code{suppress} argument be applied to?
#' @param technical list of technical arguments when models are re-estimated (see \code{\link{mirt}}
#' for details)
#' @param ... additional arguments to be passed to \code{fscores()}
#'
#' @name residuals-method
#' @aliases residuals,SingleGroupClass-method residuals,MixtureClass-method
#' residuals,MultipleGroupClass-method residuals,DiscreteClass-method
#' @rdname residuals-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#'
#' Chen, W. H. & Thissen, D. (1997). Local dependence indices for item pairs using item
#' response theory. \emph{Journal of Educational and Behavioral Statistics, 22}, 265-289.
#'
#' Edwards, M. C., Houts, C. R. & Cai, L. (2018). A Diagnostic Procedure to Detect Departures
#' From Local Independence in Item Response Theory Models.
#' \emph{Psychological Methods, 23}, 138-149.
#'
#' Yen, W. (1984). Effects of local item dependence on the fit and equating performance of the three
#' parameter logistic model. \emph{Applied Psychological Measurement, 8}, 125-145.
#' @examples
#'
#' \dontrun{
#'
#' x <- mirt(Science, 1)
#' residuals(x)
#' residuals(x, tables = TRUE)
#' residuals(x, type = 'exp')
#' residuals(x, suppress = .15)
#' residuals(x, df.p = TRUE)
#' residuals(x, df.p = TRUE, p.adjust = 'fdr') # apply FWE control
#'
#' # Pearson's X2 estimate for goodness-of-fit
#' full_table <- residuals(x, type = 'expfull')
#' head(full_table)
#' X2 <- with(full_table, sum((freq - exp)^2 / exp))
#' df <- nrow(full_table) - extract.mirt(x, 'nest') - 1
#' p <- pchisq(X2, df = df, lower.tail=FALSE)
#' data.frame(X2, df, p, row.names='Pearson-X2')
#'
#' # above FOG test as a function
#' PearsonX2 <- function(x){
#' full_table <- residuals(x, type = 'expfull')
#' X2 <- with(full_table, sum((freq - exp)^2 / exp))
#' df <- nrow(full_table) - extract.mirt(x, 'nest') - 1
#' p <- pchisq(X2, df = df, lower.tail=FALSE)
#' data.frame(X2, df, p, row.names='Pearson-X2')
#' }
#' PearsonX2(x)
#'
#'
#' # extract results manually
#' out <- residuals(x, df.p = TRUE, verbose=FALSE)
#' str(out)
#' out$df.p[1,2]
#'
#' # with and without supplied factor scores
#' Theta <- fscores(x)
#' residuals(x, type = 'Q3', Theta=Theta)
#' residuals(x, type = 'Q3', method = 'ML')
#'
#' # Edwards et al. (2018) JSI statistic
#' N <- 250
#' a <- rnorm(10, 1.7, 0.3)
#' d <- rnorm(10)
#' dat <- simdata(a, d, N=250, itemtype = '2PL')
#'
#' mod <- mirt(dat, 1)
#' residuals(mod, type = 'JSI')
#' residuals(mod, type = 'JSI', fold=FALSE) # unfolded
#'
#' # LD between items 1-2
#' aLD <- numeric(10)
#' aLD[1:2] <- rnorm(2, 2.55, 0.15)
#' a2 <- cbind(a, aLD)
#' dat <- simdata(a2, d, N=250, itemtype = '2PL')
#'
#' mod <- mirt(dat, 1)
#'
#' # JSI executed in parallel over multiple cores
#' if(interactive()) mirtCluster()
#' residuals(mod, type = 'JSI')
#'
#' }
setMethod(
f = "residuals",
signature = signature(object = 'SingleGroupClass'),
definition = function(object, type = 'LD', p.adjust = 'none',
df.p = FALSE, approx.z = FALSE,
full.scores = FALSE, QMC = FALSE,
printvalue = NULL, tables = FALSE, verbose = TRUE, Theta = NULL,
suppress = NA, theta_lim = c(-6, 6), quadpts = NULL, fold = TRUE,
upper = TRUE, technical = list(), ...)
{
dots <- list(...)
if(.hasSlot(object@Model$lrPars, 'beta'))
stop('Latent regression models not yet supported')
discrete <- use_dentype_estimate <- FALSE
if(!is.null(dots$use_dentype_estimate))
use_dentype_estimate <- dots$use_dentype_estimate
if(!is.null(dots$discrete) || use_dentype_estimate) discrete <- TRUE
K <- object@Data$K
data <- object@Data$data
N <- nrow(data)
J <- ncol(data)
nfact <- ncol(object@Fit$F)
res <- matrix(0,J,J)
diag(res) <- NA
colnames(res) <- rownames(res) <- colnames(data)
if(!is.null(Theta)){
if(!is.matrix(Theta)) Theta <- matrix(Theta)
if(nrow(Theta) > nrow(data))
Theta <- Theta[-extract.mirt(object, 'completely_missing'), , drop=FALSE]
stopifnot("Theta does not have the correct number of rows" =
nrow(Theta) == nrow(object@Data$data))
}
if(!discrete){
if(is.null(quadpts)){
if(QMC) quadpts <- 5000L
else quadpts <- object@Options$quadpts
}
if(is.nan(quadpts))
quadpts <- select_quadpts(nfact)
theta <- as.matrix(seq(theta_lim[1L], theta_lim[2L], length.out = quadpts))
if(type != 'Q3'){
Theta <- if(QMC)
QMC_quad(npts=quadpts, nfact=nfact, lim=theta_lim)
else {
if(nfact > 3L)
warning('High-dimensional models should use QMC integration instead',
call.=FALSE)
thetaComb(theta, nfact)
}
} else if(is.null(Theta)){
Theta <- fscores(object, verbose=FALSE, full.scores=TRUE, leave_missing=TRUE, ...)
}
} else {
Theta <- object@Model$Theta
}
itemnames <- colnames(data)
listtabs <- list()
calcG2 <- ifelse(type == 'LDG2', TRUE, FALSE)
if(type %in% c('LD', 'LDG2')){
if(!discrete){
groupPars <- ExtractGroupPars(object@ParObjects$pars[[object@Data$nitems + 1L]])
if(QMC){
Theta <- Theta_meanSigma_shift(Theta, groupPars$gmeans, groupPars$gcov)
prior <- rep(1/nrow(Theta), nrow(Theta))
} else {
prior <- mirt_dmvnorm(Theta,groupPars$gmeans, groupPars$gcov)
prior <- prior/sum(prior)
}
} else {
prior <- object@Internals$Prior[[1L]]
}
df <- (object@Data$K - 1) %o% (object@Data$K - 1)
diag(df) <- NA
colnames(df) <- rownames(df) <- colnames(res)
for(i in seq_len(J)){
for(j in seq_len(J)){
if(i < j){
P1 <- ProbTrace(x=object@ParObjects$pars[[i]], Theta=Theta)
P2 <- ProbTrace(x=object@ParObjects$pars[[j]], Theta=Theta)
tab <- table(data[,i], data[,j], useNA = 'no')
Etab <- matrix(0,K[i],K[j])
NN <- sum(tab)
for(k in seq_len(K[i]))
for(m in seq_len(K[j]))
Etab[k,m] <- NN * sum(P1[,k] * P2[,m] * prior)
s <- try(gamma_cor(tab) - gamma_cor(Etab), TRUE)
if(is.nan(s) || is(s, 'try-error')){
res[i,j] <- res[j,i] <- NaN
next
}
if(s == 0) s <- 1
if(calcG2){
tmp <- tab
tmp[tab == 0] <- NA
res[j,i] <- 2 * sum(tmp * log(tmp/Etab), na.rm=TRUE)
} else {
res[j,i] <- sum(((tab - Etab)^2)/Etab)
}
res[i,j] <- sign(s) * sqrt( abs(res[j,i]) / (NN * min(c(K[i],K[j]) - 1L)))
df[i,j] <- pchisq(abs(res[j,i]), df=df[j,i], lower.tail=FALSE)
if(tables){
tmp <- paste0(itemnames[i], '_', itemnames[j])
listtabs[[tmp]] <- list(Obs=tab, Exp=Etab, std_res=(tab-Etab)/sqrt(Etab))
}
if(approx.z)
res[j,i] <- sign(res[i,j]) * (sqrt(2 * res[j,i]) - sqrt(2 * df[j,i] - 1))
}
}
}
if(df.p)
df[upper.tri(df)] <- p.adjust(df[upper.tri(df)], method=p.adjust)
if(tables) return(listtabs)
if(df.p){
class(df) <- c('mirt_matrix', 'matrix')
if(verbose){
cat("Degrees of freedom (lower triangle) and p-values:\n\n")
print(df, ..., na.print = " ")
cat("\n")
}
}
if(verbose){
cat("LD matrix (lower triangle) and standardized values.\n")
cat("\nUpper triangle summary:\n")
print(round(summary(res[upper.tri(res)]), 3))
cat("\n")
}
res <- suppressMat(res, suppress=suppress, upper=upper)
class(res) <- c('mirt_matrix', 'matrix')
if(verbose) print(res, ..., na.print = " ")
if(df.p){
ret <- list(df, res)
names(ret) <- c('df.p', type)
return(invisible(ret))
}
return(invisible(res))
} else if(type == 'exp'){
r <- object@Data$Freq[[1L]]
res <- (r - object@Internals$Pl * nrow(object@Data$data)) /
sqrt(object@Internals$Pl * nrow(object@Data$data))
expected <- N * object@Internals$Pl
tabdata <- object@Data$tabdata
ISNA <- is.na(rowSums(tabdata))
expected[ISNA] <- res[ISNA] <- NA
tabdata <- data.frame(tabdata,object@Data$Freq[[1L]],expected,res)
colnames(tabdata) <- c(colnames(object@Data$tabdata),"freq","exp","std.res")
if(full.scores){
tabdata[, 'exp'] <- object@Internals$Pl / r * N
tabdata2 <- object@Data$tabdata
stabdata2 <- apply(tabdata2, 1, paste, sep='', collapse = '/')
sfulldata <- apply(object@Data$data, 1, paste, sep='', collapse = '/')
scoremat <- tabdata[match(sfulldata, stabdata2), 'exp', drop = FALSE]
res <- (1-scoremat) / sqrt(scoremat)
colnames(res) <- 'std.res'
ret <- cbind(object@Data$data, scoremat, res)
ret[is.na(rowSums(ret)), c('exp', 'std.res')] <- NA
rownames(ret) <- NULL
ret <- as.mirt_df(ret)
return(ret)
} else {
tabdata <- tabdata[do.call(order, as.data.frame(tabdata[,1:J])),]
if(!is.null(printvalue)){
if(!is.numeric(printvalue)) stop('printvalue is not a number.', call.=FALSE)
tabdata <- tabdata[abs(tabdata[ ,ncol(tabdata)]) > printvalue, ]
}
tabdata <- as.mirt_df(tabdata)
rownames(tabdata) <- NULL
return(tabdata)
}
} else if(type == 'expfull'){
K <- extract.mirt(object, 'K')
nitems <- length(K)
resp <- vector('list', nitems)
for(i in seq_len(nitems))
resp[[i]] <- 0L:(K[i]-1L)
tabdata <- expand.grid(resp)
rownames(tabdata) <- NULL
tabdata <- t(t(tabdata) + extract.mirt(object, 'mins'))
colnames(tabdata) <- colnames(extract.mirt(object, 'data'))
sv <- mod2values(object)
itemtype <- extract.mirt(object, 'itemtype')
nfact <- extract.mirt(object, 'nfact')
tmpdat <- matrix(0, nrow=2, ncol=nitems)
colnames(tmpdat) <- colnames(tabdata)
large <- mirt(tmpdat, nfact, itemtype=itemtype, pars=sv, TOL=NaN, large='return')
large$tabdata <- poly2dich(tabdata)
large$Freq$all <- rep(1L, nrow(tabdata))
large$tabdata2 <- matrix(1L)
full_object <- mirt(tabdata, nfact, itemtype=itemtype, pars=sv, TOL=NaN, large=large)
Pl <- full_object@Internals$Pl
r <- integer(length(Pl))
ro <- object@Data$Freq[[1L]]
N <- sum(ro)
for(i in seq_len(length(ro))){
pick <- colSums(t(tabdata) == object@Data$tabdata[i,]) == ncol(tabdata)
r[pick] <- ro[i]
}
res <- (r - Pl * N) / sqrt(Pl * N)
expected <- N * Pl
tabdata <- data.frame(tabdata,r,expected,res)
colnames(tabdata) <- c(colnames(object@Data$data),"freq","exp","res")
tabdata <- tabdata[do.call(order, as.data.frame(tabdata[,1:J])),]
rownames(tabdata) <- 1:length(r)
if(!is.null(printvalue)){
if(!is.numeric(printvalue)) stop('printvalue is not a number.', call.=FALSE)
tabdata <- tabdata[abs(tabdata[ ,ncol(tabdata)]) > printvalue, ]
}
tabdata <- as.mirt_df(tabdata)
return(tabdata)
} else if(type == 'Q3'){
if(discrete && !use_dentype_estimate)
stop('residual type not supported for discrete density forms', call.=FALSE)
dat <- matrix(NA, N, 2L)
diag(res) <- 1
for(i in seq_len(J)){
ei <- extract.item(object, item=i)
EI <- expected.item(ei, Theta=Theta)
dat[ ,1L] <- object@Data$data[ ,i] - EI
for(j in seq_len(J)){
if(i < j){
ej <- extract.item(object, item=j)
EJ <- expected.item(ej, Theta=Theta)
dat[,2L] <- object@Data$data[ ,j] - EJ
tmpdat <- na.omit(dat)
res[i,j] <- res[j,i] <- cor(tmpdat)[1L,2L]
}
}
}
if(verbose){
cat("Q3 summary statistics:\n")
print(round(summary(res[upper.tri(res)]), 3))
cat("\n")
}
res <- suppressMat(res, suppress=suppress, upper=upper)
class(res) <- c('mirt_matrix', 'matrix')
if(verbose) print(res, ..., na.print = " ")
return(invisible(res))
} else if(type == 'JSI'){
nfact <- extract.mirt(object, 'nfact')
stopifnot(nfact == 1L)
nitems <- extract.mirt(object, 'nitems')
technical$omp <- FALSE
as_drop <- myLapply(seq_len(nitems), function(item, mod, technical, ...){
itemtype <- extract.mirt(mod, 'itemtype')[-item]
tmpdat <- extract.mirt(mod, 'data')[,-item]
tmpmod <- mirt(tmpdat, 1L, itemtype=itemtype, SE=TRUE, verbose=FALSE,
technical=technical, ...)
ret <- sapply(coef(tmpmod, printSE=TRUE)[1:ncol(tmpdat)], function(x) x[1L:2L, 'a1'])
ret
}, mod=object, technical=technical, ...)
as <- sapply(coef(object)[1:nitems], function(x) x[1L, 'a1'])
retmat <- matrix(NA, nitems, nitems)
colnames(retmat) <- rownames(retmat) <- extract.mirt(object, 'itemnames')
for(i in seq_len(nitems)){
tmp <- as_drop[[i]]
pick <- colnames(tmp)
zs <- (as[pick] - tmp[1L, ]) / tmp[2L, ]
retmat[i, pick] <- zs
}
if(fold) retmat <- retmat + t(retmat)
retmat <- suppressMat(retmat, suppress=suppress, upper=upper)
class(retmat) <- c('mirt_matrix', 'matrix')
if(verbose){
cat("JSI summary statistics:\n")
print(round(summary(na.omit(as.vector(retmat))), 3))
cat("\n")
print(retmat, ..., na.print = " ")
}
return(invisible(retmat))
} else {
stop('specified type does not exist', call.=FALSE)
}
}
)
#' Plot various test-implied functions from models
#'
#' Plot various test implied response functions from models estimated in the mirt package.
#'
#' @param x an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{DiscreteClass}
#' @param y an arbitrary missing argument required for \code{R CMD check}
#' @param type type of plot to view. Can be
#' \describe{
#' \item{\code{'info'}}{test information function}
#' \item{\code{'rxx'}}{for the reliability function}
#' \item{\code{'infocontour'}}{for the test information contours}
#' \item{\code{'SE'}}{for the test standard error function}
#' \item{\code{'infotrace'}}{item information traceline plots}
#' \item{\code{'infoSE'}}{a combined test information and standard error plot}
#' \item{\code{'trace'}}{item probability traceline plots}
#' \item{\code{'itemscore'}}{item scoring traceline plots}
#' \item{\code{'score'}}{expected total score surface}
#' \item{\code{'scorecontour'}}{expected total score contour plot}
#' \item{\code{'posteriorTheta'}}{posterior for the latent trait distribution}
#' \item{\code{'EAPsum'}}{compares sum-scores to the expected values based
#' on the EAP for sum-scores method (see \code{\link{fscores}})}
#' }
#'
#' Note that if \code{dentype = 'empiricalhist'} was used in estimation then
#' the type \code{'empiricalhist'}
#' also will be available to generate the empirical histogram plot, and if
#' \code{dentype = 'Davidian-#'} was used then the type \code{'Davidian'}
#' will also be available to generate the curve estimates at the quadrature
#' nodes used during estimation
#' @param drop2 logical; where appropriate, for dichotomous response items drop the lowest category
#' and provide information pertaining only to the second response option?
#' @param degrees numeric value ranging from 0 to 90 used in \code{plot} to compute angle
#' for information-based plots with respect to the first dimension.
#' If a vector is used then a bubble plot is created with the summed information
#' across the angles specified (e.g., \code{degrees = seq(0, 90, by=10)})
#' @param theta_lim lower and upper limits of the latent trait (theta) to be evaluated, and is
#' used in conjunction with \code{npts}
#' @param npts number of quadrature points to be used for plotting features.
#' Larger values make plots look smoother
#' @param MI a single number indicating how many imputations to draw to form bootstrapped confidence
#' intervals for the selected test statistic. If greater than 0 a plot will be drawn with a shaded
#' region for the interval
#' @param CI a number from 0 to 1 indicating the confidence interval to select when MI input is
#' used. Default uses the 95\% confidence (CI = .95)
#' @param rot allows rotation of the 3D graphics
#' @param which.items numeric vector indicating which items to be used when plotting. Default is
#' to use all available items
#' @param facet_items logical; apply grid of plots across items? If \code{FALSE}, items will be
#' placed in one plot for each group
#' @param profile logical; provide a profile plot of response probabilities (objects returned from
#' \code{\link{mdirt}} only)
#' @param auto.key plotting argument passed to \code{\link[lattice]{lattice}}
#' @param par.strip.text plotting argument passed to \code{\link[lattice]{lattice}}
#' @param par.settings plotting argument passed to \code{\link[lattice]{lattice}}
#' @param ehist.cut a probability value indicating a threshold for excluding cases in empirical
#' histogram plots. Values larger than the default will include more points in the tails of the
#' plot, potentially squishing the 'meat' of the plot to take up less area than visually desired
#' @param main argument passed to lattice. Default generated automatically
#' @param drape logical argument passed to lattice. Default generated automatically
#' @param colorkey logical argument passed to lattice. Default generated automatically
#' @param add.ylab2 logical argument passed to lattice. Default generated automatically
#' @param ... additional arguments to be passed to lattice
#'
#' @name plot-method
#' @export
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @aliases plot,SingleGroupClass-method
#' plot,MultipleGroupClass-method plot,SingleGroupClass,missing-method
#' plot,DiscreteClass,missing-method plot,MixtureClass,missing-method
#' @rdname plot-method
#' @examples
#'
#' \dontrun{
#' x <- mirt(Science, 1, SE=TRUE)
#' plot(x)
#' plot(x, type = 'info')
#' plot(x, type = 'infotrace')
#' plot(x, type = 'infotrace', facet_items = FALSE)
#' plot(x, type = 'infoSE')
#' plot(x, type = 'rxx')
#' plot(x, type = 'posteriorTheta')
#'
#' # confidence interval plots when information matrix computed
#' plot(x)
#' plot(x, MI=100)
#' plot(x, type='info', MI=100)
#' plot(x, type='SE', MI=100)
#' plot(x, type='rxx', MI=100)
#'
#' # use the directlabels package to put labels on tracelines
#' library(directlabels)
#' plt <- plot(x, type = 'trace')
#' direct.label(plt, 'top.points')
#'
#' # additional modifications can be made via update().
#' # See ?update.trellis for further documentation
#' plt
#' update(plt, ylab = expression(Prob(theta)),
#' main = "Item Traceline Functions") # ylab/main changed
#'
#' set.seed(1234)
#' group <- sample(c('g1','g2'), nrow(Science), TRUE)
#' x2 <- multipleGroup(Science, 1, group)
#' plot(x2)
#' plot(x2, type = 'trace')
#' plot(x2, type = 'trace', which.items = 1:2)
#' plot(x2, type = 'itemscore', which.items = 1:2)
#' plot(x2, type = 'trace', which.items = 1, facet_items = FALSE) #facet by group
#' plot(x2, type = 'info')
#'
#' x3 <- mirt(Science, 2)
#' plot(x3, type = 'info')
#' plot(x3, type = 'SE', theta_lim = c(-3,3))
#'
#' }
setMethod(
f = "plot",
signature = signature(x = 'SingleGroupClass', y = 'missing'),
definition = function(x, y, type = 'score', npts = 200, drop2 = TRUE, degrees = 45,
theta_lim = c(-6,6), which.items = 1:extract.mirt(x, 'nitems'),
MI = 0, CI = .95, rot = list(xaxis = -70, yaxis = 30, zaxis = 10),
facet_items = TRUE, main = NULL,
drape = TRUE, colorkey = TRUE, ehist.cut = 1e-10, add.ylab2 = TRUE,
par.strip.text = list(cex = 0.7),
par.settings = list(strip.background = list(col = '#9ECAE1'),
strip.border = list(col = "black")),
auto.key = list(space = 'right', points=FALSE, lines=TRUE),
profile = FALSE, ...)
{
dots <- list(...)
if(!(type %in% c('info', 'SE', 'infoSE', 'rxx', 'trace', 'score', 'itemscore',
'infocontour', 'infotrace', 'scorecontour', 'empiricalhist', 'Davidian',
'EAPsum', 'posteriorTheta')))
stop('type supplied is not supported')
if (any(degrees > 90 | degrees < 0))
stop('Improper angle specified. Must be between 0 and 90.', call.=FALSE)
rot <- list(x = rot[[1]], y = rot[[2]], z = rot[[3]])
nfact <- x@Model$nfact
if(nfact > 3) stop("Can't plot high dimensional solutions.", call.=FALSE)
J <- x@Data$nitems
if(x@ParObjects$pars[[J + 1L]]@dentype == 'custom')
theta_lim <- x@Internals$theta_lim
theta <- seq(theta_lim[1L],theta_lim[2L],length.out=npts/(nfact^2))
ThetaFull <- Theta <- thetaComb(theta, nfact)
prodlist <- attr(x@ParObjects$pars, 'prodlist')
if(all(x@Data$K[which.items] == 2L) && facet_items) auto.key <- FALSE
if(length(prodlist))
ThetaFull <- prodterms(Theta,prodlist)
if(length(degrees) > ncol(ThetaFull)) type <- 'infoangle'
if(length(degrees) == 1L) degrees <- rep(degrees, ncol(ThetaFull))
info <- numeric(nrow(ThetaFull))
if(type %in% c('info', 'infocontour', 'rxx', 'SE', 'infoSE'))
info <- testinfo(x, ThetaFull, degrees = degrees, which.items=which.items)
if(type == 'infoangle'){
for(i in seq_len(length(degrees)))
info <- info + testinfo(x, ThetaFull, degrees = rep(degrees[i], ncol(ThetaFull)),
which.items=which.items)
}
mins <- x@Data$mins
maxs <- extract.mirt(x, 'K') + mins - 1
rotate <- if(is.null(dots$rotate)) 'none' else dots$rotate
if (x@Options$exploratory){
if(!is.null(dots$rotate)){
so <- summary(x, verbose=FALSE, digits=5, ...)
a <- rotateLambdas(so) * 1.702
for(i in 1:J)
x@ParObjects$pars[[i]]@par[1:nfact] <- a[i, ]
}
}
itemtrace <- computeItemtrace(x@ParObjects$pars, ThetaFull, x@Model$itemloc,
CUSTOM.IND=x@Internals$CUSTOM.IND)
score <- c()
for(i in 1:J)
score <- c(score, (0:(x@Data$K[i]-1) + mins[i]) * (i %in% which.items))
score <- matrix(score, nrow(itemtrace), ncol(itemtrace), byrow = TRUE)
plt <- data.frame(cbind(info,score=rowSums(score*itemtrace),Theta=Theta))
bundle <- length(which.items) != J
gp <- ExtractGroupPars(x@ParObjects$pars[[J+1]])
if(MI > 0L && nfact == 1L){
tmpx <- x
if(!x@Options$SE)
stop('Must compute an information matrix', call.=FALSE)
covB <- x@vcov
pre.ev <- eigen(covB)
names <- colnames(covB)
tmp <- lapply(names, function(x, split){
as.numeric(strsplit(x, split=split)[[1L]][-1L])
}, split='\\.')
imputenums <- do.call(c, tmp)
CIscore <- CIinfo <- CIrxx <- matrix(0, MI, length(plt$score))
for(i in seq_len(MI)){
while(TRUE){
tmp <- try(imputePars(pars=x@ParObjects$pars, pre.ev=pre.ev,
imputenums=imputenums, constrain=x@Model$constrain),
silent=TRUE)
if(!is(tmp, 'try-error')) break
}
tmpx@ParObjects$pars <- tmp
gp2 <- ExtractGroupPars(tmp[[J+1]])
itemtrace <- computeItemtrace(tmpx@ParObjects$pars, ThetaFull, x@Model$itemloc,
CUSTOM.IND=x@Internals$CUSTOM.IND)
tmpscore <- rowSums(score * itemtrace)
CIscore[i, ] <- tmpscore
CIinfo[i, ] <- testinfo(tmpx, ThetaFull)
CIrxx[i, ] <- CIinfo[i, ] / (CIinfo[i, ] + 1/gp2$gcov[1L,1L])
}
}
mins <- mins[which.items]
maxs <- maxs[which.items]
ybump <- (max(maxs) - min(mins))/15
ybump_full <- (sum(maxs) - sum(mins))/15
if(type == 'EAPsum'){
if(is.null(main))
main <- "Expected vs Observed Sum-Scores"
fs <- fscores(x, method = 'EAPsum', full.scores=FALSE, verbose=FALSE, ...)
plt <- with(fs, data.frame(Scores=Sum.Scores, y=c(observed, expected),
group = rep(c('observed', 'expected'), each=nrow(fs))))
return(xyplot(y~Scores, plt, type='l', main = main, group=plt$group,
auto.key=auto.key, xlab = expression(Sum-Score), ylab=expression(n),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
if(nfact == 3){
colnames(plt) <- c("info", "score", "Theta1", "Theta2", "Theta3")
plt$SE <- 1 / sqrt(plt$info)
if(type == 'infocontour'){
if(is.null(main)){
main <- paste("Test Information Contour")
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(contourplot(info ~ Theta1 * Theta2 | Theta3, data = plt,
main = main, xlab = expression(theta[1]),
ylab = expression(theta[2]),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'scorecontour'){
if(is.null(main)){
main <- paste("Expected Score Contour")
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(contourplot(score ~ Theta1 * Theta2 | Theta3, data = plt,
ylim=c(sum(mins)-ybump_full, sum(maxs)+ybump_full),
main = main, xlab = expression(theta[1]),
ylab = expression(theta[2]), ylim=c(sum(mins)-.1, sum(maxs)+.1),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'info'){
if(is.null(main)){
main <- "Test Information"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(wireframe(info ~ Theta1 + Theta2 | Theta3, data = plt, main = main,
zlab=expression(I(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'SEcontour'){
if(is.null(main)){
main <- "Test Standard Errors"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(contourplot(score ~ Theta1 * Theta2 | Theta3, data = plt,
main = main, xlab = expression(theta[1]),
ylab = expression(theta[2]),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'score'){
if(is.null(main)){
main <- if(bundle) "Expected Bundle Score" else "Expected Total Score"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(wireframe(score ~ Theta1 + Theta2 | Theta3, data = plt, main = main,
ylim=c(sum(mins)-ybump_full, sum(maxs)+ybump_full),
zlab=expression(Total(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings,
ylim=c(sum(mins)-.1, sum(maxs)+.1), ...))
} else if(type == 'SE'){
if(is.null(main)){
main <- "Test Standard Errors"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(wireframe(SE ~ Theta1 + Theta2 | Theta3, data = plt, main = main,
zlab=expression(SE(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
stop('plot type not supported for three dimensional model', call.=FALSE)
}
} else if(nfact == 2){
colnames(plt) <- c("info", "score", "Theta1", "Theta2")
plt$SE <- 1 / sqrt(plt$info)
if(type == 'infocontour'){
if(is.null(main)){
main <- paste("Test Information Contour")
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(contourplot(info ~ Theta1 * Theta2, data = plt,
main = main, xlab = expression(theta[1]),
ylab = expression(theta[2]),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'scorecontour'){
if(is.null(main)){
main <- paste("Expected Score Contour")
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(contourplot(score ~ Theta1 * Theta2, data = plt,
ylim=c(sum(mins)-ybump_full, sum(maxs)+ybump_full),
main = main, xlab = expression(theta[1]),
ylab = expression(theta[2]),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'info'){
if(is.null(main)){
main <- "Test Information"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(wireframe(info ~ Theta1 + Theta2, data = plt, main = main,
zlab=expression(I(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'trace'){
if(is.null(main))
main <- 'Item Probability Functions'
P <- vector('list', length(which.items))
names(P) <- colnames(x@Data$data)[which.items]
ind <- 1L
for(i in which.items){
tmp <- probtrace(extract.item(x, i), ThetaFull)
if(ncol(tmp) == 2L && drop2) tmp <- tmp[,2, drop=FALSE]
tmp2 <- data.frame(P=as.numeric(tmp), cat=gl(ncol(tmp), k=nrow(ThetaFull),
labels=paste0('P', seq_len(ncol(tmp)))))
P[[ind]] <- tmp2
ind <- ind + 1L
}
nrs <- sapply(P, nrow)
Pstack <- do.call(rbind, P)
names <- c()
for(i in seq_len(length(nrs)))
names <- c(names, rep(names(P)[i], nrs[i]))
plotobj <- data.frame(Pstack, item=names, Theta=ThetaFull)
plotobj$item <- factor(plotobj$item, levels = colnames(x@Data$data)[which.items])
return(wireframe(P ~ Theta.1 * Theta.2 | item, plotobj, zlim = c(-0.1,1.1), group=plotobj$cat,
zlab=expression(P(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'infotrace'){
if(is.null(main))
main <- 'Item Information'
I <- matrix(NA, nrow(Theta), J)
for(i in which.items)
I[,i] <- iteminfo(extract.item(x, i), ThetaFull, degrees=degrees)
I <- t(na.omit(t(I)))
items <- rep(colnames(x@Data$data)[which.items], each=nrow(Theta))
plotobj <- data.frame(I = as.numeric(I), Theta=ThetaFull, item=items)
plotobj$item <- factor(plotobj$item, levels = colnames(x@Data$data)[which.items])
return(wireframe(I ~ Theta.1 * Theta.2 | item, plotobj, group=plotobj$cat,
zlab=expression(P(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'itemscore'){
if(is.null(main))
main <- 'Expected Item Score'
S <- vector('list', length(which.items))
names(S) <- colnames(x@Data$data)[which.items]
ind <- 1L
for(i in which.items){
S[[ind]] <- expected.item(extract.item(x, i), ThetaFull, mins[i])
ind <- ind + 1L
}
Sstack <- do.call(c, S)
names <- rep(names(S), each = nrow(ThetaFull))
plotobj <- data.frame(S=Sstack, item=names, Theta=ThetaFull)
plotobj$item <- factor(plotobj$item, levels = colnames(x@Data$data)[which.items])
return(wireframe(S ~ Theta.1 * Theta.2 | item, plotobj, group=plotobj$cat,
zlab=expression(P(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'SEcontour'){
if(is.null(main)){
main <- "Test Standard Errors"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(contourplot(score ~ Theta1 * Theta2, data = plt,
main = main, xlab = expression(theta[1]),
ylab = expression(theta[2]),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'score'){
if(is.null(main)){
main <- if(bundle) "Expected Bundle Score" else "Expected Total Score"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(wireframe(score ~ Theta1 + Theta2, data = plt, main = main,
zlim=c(sum(mins)-ybump_full, sum(maxs)+ybump_full),
zlab=expression(Total(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'infoangle'){
if(is.null(main)){
main <- 'Information across different angles'
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
graphics::symbols(plt$Theta1, plt$Theta2, circles = sqrt(plt$info/pi), inches = .35, fg='white', bg='blue',
xlab = expression(theta[1]), ylab = expression(theta[2]),
main = main)
} else if(type == 'SE'){
if(is.null(main)){
main <- "Test Standard Errors"
if(x@Options$exploratory) main <- paste0(main, ' (rotate = \'', rotate, '\')')
}
return(wireframe(SE ~ Theta1 + Theta2, data = plt, main = main,
zlab=expression(SE(theta)), xlab=expression(theta[1]), ylab=expression(theta[2]),
scales = list(arrows = FALSE), screen = rot, colorkey = colorkey, drape = drape,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
stop('plot type not supported for two dimensional model', call.=FALSE)
}
} else {
colnames(plt) <- c("info", "score", "Theta")
plt$SE <- 1 / sqrt(plt$info)
plt$rxx <- plt$info / (plt$info + 1/gp$gcov[1L,1L])
if(MI > 0){
bs_range <- function(x, CI){
ss <- sort(x)
N <- length(ss)
ret <- c(upper = ss[ceiling(N * (1 - (1-CI)/2))],
middle = median(x),
lower = ss[floor(N * (1-CI)/2)])
ret
}
tmp <- apply(CIscore, 2, bs_range, CI=CI)
plt$CIscoreupper <- tmp['upper', ]
plt$CIscorelower <- tmp['lower', ]
tmp <- apply(CIinfo, 2, bs_range, CI=CI)
plt$CIinfoupper <- tmp['upper', ]
plt$CIinfolower <- tmp['lower', ]
plt$CISElower <- 1/sqrt(tmp['upper', ])
plt$CISEupper <- 1/sqrt(tmp['lower', ])
tmp <- apply(CIrxx, 2, bs_range, CI=CI)
plt$CIrxxupper <- tmp['upper', ]
plt$CIrxxlower <- tmp['lower', ]
}
if(type == 'info'){
if(is.null(main))
main <- 'Test Information'
if(MI > 0){
return(xyplot(info ~ Theta, data=plt,
upper=plt$CIinfoupper, lower=plt$CIinfolower,
panel = function(x, y, lower, upper, ...){
panel.polygon(c(x, rev(x)), c(upper, rev(lower)),
col="#E6E6E6", border = FALSE, ...)
panel.xyplot(x, y, type='l', lty=1,...)
},
main = main, ylim=c(min(plt$CIinfolower), max(plt$CIinfoupper)),
ylab = expression(I(theta)), xlab = expression(theta),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
return(xyplot(info~Theta, plt, type='l', main = main,
xlab = expression(theta), ylab=expression(I(theta)),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
} else if(type == 'rxx'){
if(is.null(main))
main <- 'Reliability'
if(MI > 0){
return(xyplot(rxx ~ Theta, data=plt,
upper=plt$CIrxxupper, lower=plt$CIrxxlower,
panel = function(x, y, lower, upper, ...){
panel.polygon(c(x, rev(x)), c(upper, rev(lower)),
col="#E6E6E6", border = FALSE, ...)
panel.xyplot(x, y, type='l', lty=1,...)
},
main = main, ylim=c(-0.1, 1.1),
ylab = expression(r[xx](theta)), xlab = expression(theta),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
return(xyplot(rxx~Theta, plt, type='l', main = main, ylim=c(-0.1, 1.1),
xlab = expression(theta), ylab=expression(r[xx](theta)),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
} else if(type == 'SE'){
if(is.null(main))
main <- 'Test Standard Errors'
if(MI > 0){
return(xyplot(SE ~ Theta, data=plt,
upper=plt$CISEupper, lower=plt$CISElower,
panel = function(x, y, lower, upper, ...){
panel.polygon(c(x, rev(x)), c(upper, rev(lower)),
col="#E6E6E6", border = FALSE, ...)
panel.xyplot(x, y, type='l', lty=1,...)
},
main = main, ylim=c(min(plt$CISElower), max(plt$CISEupper)),
ylab = expression(I(theta)), xlab = expression(theta),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
return(xyplot(SE~Theta, plt, type='l', main = main,
xlab = expression(theta), ylab=expression(SE(theta)),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
} else if(type == 'infoSE'){
if(is.null(main))
main <- 'Test Information and Standard Errors'
par.settings <- c(par.settings,
lattice::simpleTheme(col = c("#0080ff",'red'), lty = 1:2))
obj1 <- xyplot(info~Theta, plt, type='l', main = main,
xlab = expression(theta), ylab=expression(I(theta)),
par.strip.text=par.strip.text, par.settings=par.settings)
obj2 <- xyplot(SE~Theta, plt, type='l', ylab=expression(SE(theta)),
par.strip.text=par.strip.text, par.settings=par.settings)
if(requireNamespace("latticeExtra", quietly = TRUE)){
return(latticeExtra::doubleYScale(obj1, obj2, add.ylab2 = add.ylab2, ...))
} else {
stop('latticeExtra package is not available. Please install.', call.=FALSE)
}
} else if(type == 'trace'){
if(is.null(main))
main <- 'Item Probability Functions'
P <- vector('list', length(which.items))
names(P) <- colnames(x@Data$data)[which.items]
ind <- 1L
alltwocats <- all(extract.mirt(x, 'K')[which.items] == 2L)
for(i in which.items){
tmp <- probtrace(extract.item(x, i), ThetaFull)
if(ncol(tmp) == 2L && (facet_items || (!facet_items && alltwocats)) && drop2)
tmp <- tmp[,2, drop=FALSE]
tmp2 <- data.frame(P=as.numeric(tmp), cat=gl(ncol(tmp), k=nrow(Theta),
labels=paste0('P', seq_len(ncol(tmp)))))
P[[ind]] <- tmp2
ind <- ind + 1L
}
nrs <- sapply(P, nrow)
Pstack <- do.call(rbind, P)
names <- c()
for(i in seq_len(length(nrs)))
names <- c(names, rep(names(P)[i], nrs[i]))
plotobj <- data.frame(Pstack, item=names, Theta=Theta)
plotobj$item <- factor(plotobj$item, levels = colnames(x@Data$data)[which.items])
if(facet_items){
return(xyplot(P ~ Theta|item, plotobj, ylim = c(-0.1,1.1), groups = cat,
xlab = expression(theta), ylab = expression(P(theta)),
auto.key = auto.key, type = 'l', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
return(xyplot(P ~ Theta|cat, plotobj, groups=plotobj$item, ylim = c(-0.1,1.1),
xlab = expression(theta), ylab = expression(P(theta)),
auto.key = auto.key, type = 'l', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
} else if(type == 'itemscore'){
if(is.null(main))
main <- 'Expected Item Score'
S <- vector('list', length(which.items))
names(S) <- colnames(x@Data$data)[which.items]
ind <- 1L
for(i in which.items){
S[[ind]] <- expected.item(extract.item(x, i), ThetaFull, mins[i])
ind <- ind + 1L
}
Sstack <- do.call(c, S)
names <- rep(names(S), each = nrow(ThetaFull))
plotobj <- data.frame(S=Sstack, item=names, Theta=Theta)
plotobj$item <- factor(plotobj$item, levels = colnames(x@Data$data)[which.items])
if(facet_items){
return(xyplot(S ~ Theta|item, plotobj, ylim=c(min(mins)-ybump, max(maxs)+ybump),
xlab = expression(theta), ylab = expression(S(theta)),
auto.key = auto.key, type = 'l', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
return(xyplot(S ~ Theta, plotobj, groups=plotobj$item, ylim=c(min(mins)-.1, max(maxs)+.1),
xlab = expression(theta), ylab = expression(S(theta)),
auto.key = auto.key, type = 'l', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
} else if(type == 'infotrace'){
if(is.null(main))
main <- 'Item Information'
I <- matrix(NA, nrow(Theta), J)
for(i in which.items)
I[,i] <- iteminfo(extract.item(x, i), ThetaFull)
I <- t(na.omit(t(I)))
items <- rep(colnames(x@Data$data)[which.items], each=nrow(Theta))
plotobj <- data.frame(I = as.numeric(I), Theta=Theta, item=items)
plotobj$item <- factor(plotobj$item, levels = colnames(x@Data$data)[which.items])
if(facet_items){
return(xyplot(I ~ Theta|item, plotobj,
xlab = expression(theta), ylab = expression(I(theta)),
auto.key = auto.key, type = 'l', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
return(xyplot(I ~ Theta, plotobj, groups = plotobj$item,
xlab = expression(theta), ylab = expression(I(theta)),
auto.key = auto.key, type = 'l', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
} else if(type == 'score'){
if(is.null(main))
main <- if(bundle) "Expected Bundle Score" else "Expected Total Score"
if(MI > 0){
return(xyplot(score ~ Theta, data=plt,
ylim=c(sum(mins)-ybump_full, sum(maxs)+ybump_full),
upper=plt$CIscoreupper, lower=plt$CIscorelower,
panel = function(x, y, lower, upper, ...){
panel.polygon(c(x, rev(x)), c(upper, rev(lower)),
col="#E6E6E6", border = FALSE, ...)
panel.xyplot(x, y, type='l', lty=1,...)
},
main = main,
ylab = expression(T(theta)), xlab = expression(theta),
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else {
return(xyplot(score ~ Theta, plt, ylim=c(sum(mins)-ybump_full, sum(maxs)+ybump_full),
xlab = expression(theta), ylab = expression(T(theta)),
type = 'l', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}
} else if(type == 'empiricalhist'){
if(!(x@Options$dentype %in% c('EH', 'EHW')))
stop('Empirical histogram was not estimated for this object', call.=FALSE)
if(is.null(main))
main <- 'Empirical Histogram'
Prior <- x@Internals$Prior[[1L]]
Theta <- x@Model$Theta
# Prior <- Prior * nrow(x@Data$data)
cuts <- cut(Theta, floor(npts/2))
Prior <- do.call(c, lapply(split(Prior, cuts), mean))
Theta <- do.call(c, lapply(split(Theta, cuts), mean))
keep1 <- min(which(Prior > ehist.cut))
keep2 <- max(which(Prior > ehist.cut))
plt <- data.frame(Theta = Theta, Prior = Prior)
plt <- plt[keep1:keep2, , drop=FALSE]
return(xyplot(Prior ~ Theta, plt,
xlab = expression(theta), ylab = 'Density',
type = 'b', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'Davidian'){
if(x@Options$dentype != 'Davidian')
stop('Davidian curve was not estimated for this object', call.=FALSE)
if(is.null(main))
main <- 'Davidian Curve'
Prior <- x@Internals$Prior[[1L]]
Theta <- x@Model$Theta
# Prior <- Prior * nrow(x@Data$data)
plt <- data.frame(Theta = Theta, Prior = Prior)
return(xyplot(Prior ~ Theta, plt,
xlab = expression(theta), ylab = 'Density',
type = 'b', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
} else if(type == 'posteriorTheta'){
if(is.null(main))
main <- 'Posterior Distribution of Latent Trait'
plt <- extract.mirt(x, 'thetaPosterior')[[1]]
return(xyplot(posterior ~ Theta, plt,
xlab = expression(theta), ylab = 'Density',
type = 'b', main = main,
par.strip.text=par.strip.text, par.settings=par.settings, ...))
}else {
stop('plot not supported for unidimensional models', call.=FALSE)
}
}
}
)
mirt2traditional <- function(x, vcov, nfact){
cls <- class(x)
opar <- par <- x@par
which.a <- which(x@par[1L:nfact] != 0)
if(length(which.a) != 1L) return(x)
a.nms <- if(nfact == 1L) 'a' else paste0('a', 1L:nfact)
if(cls != 'GroupPars')
ncat <- x@ncat
if(cls == 'dich'){
fns <- vector('list', nfact + 3L)
fns[[nfact+1L]] <- function(par, index, opar){
if(index == (nfact + 1L)){
opar[c(which.a, nfact + 1L)] <- par
ret <- -opar[nfact + 1L]/opar[which.a]
}
ret
}
fns[[nfact+2L]] <- function(par, index, opar){
if(index == nfact + 2L)
ret <- plogis(par)
ret
}
fns[[nfact+3L]] <- function(par, index, opar){
if(index == nfact + 3L)
ret <- plogis(par)
ret
}
delta_index <- c(as.list(rep(NA, nfact)),
list(c(which.a, nfact + 1L),
nfact + 2L, nfact+3L))
par[nfact + 1L] <- -par[nfact + 1L]/par[which.a]
par[nfact + 2L] <- plogis(par[nfact + 2L])
par[nfact + 3L] <- plogis(par[nfact + 3L])
names(par) <- c(a.nms, 'b', 'g', 'u')
} else if(cls == 'fivePL'){
fns <- vector('list', nfact + 4L)
fns[[nfact+1L]] <- function(par, index, opar){
if(index == (nfact + 1L)){
opar[c(which.a, nfact + 1L)] <- par
ret <- -opar[nfact + 1L]/opar[which.a]
}
ret
}
fns[[nfact+2L]] <- function(par, index, opar){
if(index == nfact + 2L)
ret <- plogis(par)
ret
}
fns[[nfact+3L]] <- function(par, index, opar){
if(index == nfact + 3L)
ret <- plogis(par)
ret
}
fns[[nfact+4L]] <- function(par, index, opar){
if(index == nfact + 4L)
ret <- exp(par)
ret
}
delta_index <- c(as.list(rep(NA, nfact)),
list(c(which.a, nfact + 1L),
nfact + 2L, nfact+3L, nfact+4L))
par[nfact + 1L] <- -par[nfact + 1L]/par[which.a]
par[nfact + 2L] <- plogis(par[nfact + 2L])
par[nfact + 3L] <- plogis(par[nfact + 3L])
par[nfact + 4L] <- exp(par[nfact + 4L])
names(par) <- c(a.nms, 'b', 'g', 'u', 'S')
} else if(cls == 'graded'){
fns <- vector('list', ncat + nfact-1L)
for(i in 2L:ncat - 1L){
fns[[i + nfact]] <- function(par, index, opar){
if(index > nfact){
opar[c(which.a, index)] <- par
ret <- -opar[index]/opar[which.a]
}
ret
}
}
delta_index <- vector('list', ncat + nfact - 1L)
for(i in 1:nfact)
delta_index[[i]] <- NA
for(i in 2:ncat-1L){
par[i+nfact] <- -par[i+nfact]/par[which.a]
delta_index[[i+nfact]] <- c(which.a, i+nfact)
}
names(par) <- c(a.nms, paste0('b', 2:ncat-1L))
} else if(cls == 'gpcm'){
fns <- vector('list', ncat+nfact)
for(i in 2L:ncat-1L){
fns[[i+nfact]] <- function(par, index, opar){
if(index > nfact){
if(index == (nfact+1))
opar[c(which.a, ncat + nfact + 2L)] <- par
else opar[c(which.a, ncat + index + nfact - 1L,
ncat + index + nfact)] <- par
par <- opar
ds <- par[(nfact+1):length(par)]/par[which.a]
ds <- ds[-seq_len(ncat)]
newd <- numeric(length(ds)-1L)
for(i in 2:length(ds))
newd[i-1L] <- -(ds[i] - ds[i-1L])
ret <- c(par[1:nfact], newd)
ret <- ret[index]
}
ret
}
}
delta_index <- vector('list', ncat+nfact-1L)
for(i in 1:nfact)
delta_index[[i]] <- NA
ds <- par[(nfact+1):length(par)]/par[which.a]
ds <- ds[-seq_len(ncat)]
newd <- numeric(length(ds)-1L)
for(i in 2:length(ds))
newd[i-1L] <- -(ds[i] - ds[i-1L])
tmp <- rbind(1:(ncat-1), 2:ncat)
for(i in 2:length(ds)){
newd[i-1L] <- -(ds[i] - ds[i-1L])
delta_index[[i+nfact-1L]] <- c(which.a, tmp[,i-1L] + ncat + nfact)
}
delta_index[[nfact+1L]] <- c(which.a, ncat + nfact + 2L)
par <- c(x@par[1:nfact], newd)
names(par) <- c(a.nms, paste0('b', 1:length(newd)))
x@est <- x@est[c(1:nfact, (ncat+nfact+2L):length(x@est))]
} else if(cls == 'nominal'){
if(nfact > 1L) return(x)
fns <- vector('list', ncat*2)
for(i in 2L:length(par)-1L){
fns[[i]] <- function(par, index, opar){
if(index <= floor(length(opar)/2)){
as <- par[2:(ncat+1)] * par[1]
as <- as - mean(as)
ret <- as[index]
} else {
ds <- par
ds <- ds - mean(ds)
ret <- ds[index-ncat]
}
ret
}
}
delta_index <- vector('list', ncat*2)
for(i in 1L:(ncat)) delta_index[[i]] <- 1L:(ncat+1L)
for(i in (ncat+1L):(ncat*2)) delta_index[[i]] <- (ncat+2L):length(par)
as <- par[2:(ncat+1)] * par[1]
as <- as - mean(as)
ds <- par[(ncat+2):length(par)]
ds <- ds - mean(ds)
par <- c(as, ds)
names(par) <- c(paste0('a', 1:ncat), paste0('c', 1:ncat))
x@est <- rep(TRUE, ncat*2)
x@SEpar <- rep(as.numeric(NA), ncat*2)
} else if(cls == 'nestlogit'){
if(nfact > 1L) return(x)
fns <- vector('list', ncat*2 + 4)
fns[[2]] <- function(par, index, opar){
if(index == 2L){
opar[1L:2L] <- par
ret <- -opar[2L]/opar[1L]
}
ret
}
fns[[3]] <- function(par, index, opar){
if(index == 3L)
ret <- plogis(par)
ret
}
fns[[4]] <- function(par, index, opar){
if(index == 4L)
ret <- plogis(par)
ret
}
for(i in (2L:length(par)-1L) + 4){
fns[[i]] <- function(par, index, opar){
if(index <= (5 + ncat-2)){
as <- par
as <- as - mean(as)
ret <- as[index-4]
} else {
ds <- par
ds <- ds - mean(ds)
ret <- ds[index-4-(ncat-1)]
}
ret
}
}
delta_index <- vector('list', (ncat-1)*2)
for(i in 1L:(ncat-1)) delta_index[[i]] <- 1L:(ncat-1) + 4
for(i in ncat:((ncat-1)*2)) delta_index[[i]] <- 1L:(ncat-1) + 4 + ncat - 1L
delta_index <- c(list(NA, 1L:2L, 3L, 4L), delta_index)
par1 <- par[1:4]
par1[2] <- -par1[2]/par1[1]
par1[3] <- plogis(par1[3])
par1[4] <- plogis(par1[4])
names(par1) <- c(a.nms, 'b', 'g', 'u')
par2 <- par[5:length(par)]
as <- par2[1:(ncat-1)]
as <- as - mean(as)
ds <- par2[-c(1:(ncat-1))]
ds <- ds - mean(ds)
names(as) <- paste0('a', 1:(ncat-1))
names(ds) <- paste0('c', 1:(ncat-1))
par <- c(par1, as, ds)
x@est <- c(x@est[1:4], rep(TRUE, (ncat-1)*2))
x@SEpar <- c(x@SEpar[1], as.numeric(rep(NA, (ncat-1)*2 + 3)))
}
x@par <- par
names(x@est) <- names(par)
x@parnames <- names(x@par)
if(length(vcov) == 0L || (is.na(vcov[1L,1L]) || !(cls %in%
c('dich', 'graded', 'gpcm', 'nominal', 'nestlogit')))){
x@SEpar <- numeric()
} else {
nms <- colnames(vcov)
splt <- strsplit(nms, "\\.")
splt <- lapply(splt, function(x) as.integer(x[-1L]))
for(i in 1L:length(delta_index)){
if(!x@est[i]) next
if(is.na(delta_index[[i]][1L])) next
grad <- numerical_deriv(opar[delta_index[[i]]], fns[[i]], opar=opar, index=i)
parnum <- x@parnum[delta_index[[i]]]
pick <- numeric(length(grad))
for(j in 1L:length(parnum))
pick[j] <- suppressWarnings(min(which(do.call(c, lapply(splt, function(x)
any(x %in% parnum[j]))))))
grad <- grad[is.finite(pick)]
pick <- pick[is.finite(pick)]
x@SEpar[i] <- as.vector(sqrt(grad %*% vcov[pick, pick, drop=FALSE] %*% grad))
}
if(cls == 'gpcm') x@SEpar <- x@SEpar[1L:length(x@par)]
}
x
}
#' Convert traditional IRT metric into slope-intercept form used in mirt
#'
#' This is a helper function for users who have previously available traditional/classical
#' IRT parameters and want to know the equivalent slope-intercept translation used in \code{mirt}.
#' Note that this function assumes that the supplied models are unidimensional by definition (i.e.,
#' will have only one slope/discrimination) and in the logistic metric (i.e., logistic-ogive
#' scaling coefficient D=1). If there is no supported slope-intercept transformation
#' available then the original vector of parameters will be returned by default.
#'
#' Supported class transformations for the \code{cls} input are:
#'
#' \describe{
#' \item{Rasch, 2PL, 3PL, 3PLu, 4PL}{
#' Form must be: (discrimination, difficulty, lower-bound, upper-bound)
#' }
#' \item{graded}{
#' Form must be: (discrimination, difficulty 1, difficulty 2, ..., difficulty k-1)
#' }
#' \item{gpcm}{
#' Form must be: (discrimination, difficulty 1, difficulty 2, ..., difficulty k-1)
#' }
#' \item{nominal}{
#' Form must be: (discrimination 1, discrimination 2, ..., discrimination k,
#' difficulty 1, difficulty 2, ..., difficulty k)
#' }
#' }
#'
#' @param x a vector of parameters to transform
#' @param cls the class or itemtype of the supplied model
#' @param ncat the number of categories implied by the IRT model
#'
#' @return a named vector of slope-intercept parameters (if supported)
#' @export
#'
#' @examples
#'
#' # classical 3PL model
#' vec <- c(a=1.5, b=-1, g=.1, u=1)
#' slopeint <- traditional2mirt(vec, '3PL', ncat=2)
#' slopeint
#'
#' # classical graded model (four category)
#' vec <- c(a=1.5, b1=-1, b2=0, b3=1.5)
#' slopeint <- traditional2mirt(vec, 'graded', ncat=4)
#' slopeint
#'
#' # classical generalize partial credit model (four category)
#' vec <- c(a=1.5, b1=-1, b2=0, b3=1.5)
#' slopeint <- traditional2mirt(vec, 'gpcm', ncat=4)
#' slopeint
#'
#' # classical nominal model (4 category)
#' vec <- c(a1=.5, a2 = -1, a3=1, a4=-.5, d1=1, d2=-1, d3=-.5, d4=.5)
#' slopeint <- traditional2mirt(vec, 'nominal', ncat=4)
#' slopeint
#'
#'
traditional2mirt <- function(x, cls, ncat){
cls <- toInternalItemtype(cls)
if(cls == 'dich'){
par <- x
par[2L] <- -par[2L]*par[1L]
names(par) <- c('a1', 'd', 'g', 'u')
} else if(cls == 'graded'){
par <- x
for(i in 2L:ncat)
par[i] <- -par[i]*par[1L]
names(par) <- c('a1', paste0('d', 1:(length(par)-1)))
} else if(cls %in% c('gpcm', 'gpcmIRT')){
if(cls == 'gpcmIRT'){
x[-c(1, length(x))] <- x[-c(1, length(x))] - x[length(x)]
x <- x[-length(x)]
}
par <- c(x[1L], 0L:(ncat-1L), 0, x[-1L])
ds <- -par[-c(1:(ncat+1))]*par[1]
newd <- numeric(length(ds))
for(i in length(ds):2L)
newd[i] <- (ds[i] + ds[i-1L])
for(i in length(newd):3L){
pick <- if(i %% 2 == 0) seq(2, i, by = 2) else seq(1, i, by = 2)
newd[i] <- sum(newd[pick])
}
par <- c(par[1:(ncat+1)], newd)
names(par) <- c('a1', paste0('ak', 0:(ncat-1)), paste0('d', 0:(ncat-1)))
} else if(cls == 'nominal'){
as <- x[seq_len(length(x)/2)]
ds <- x[-seq_len(length(x)/2)]
a1 <- (as[ncat] - as[1L]) / (ncat-1L)
ak <- 1:ncat - 1
for(i in 2:(ncat-1))
ak[i] <- -(as[1L] - as[i]) / a1
dk <- ak
for(i in 2:ncat)
dk[i] <- ds[i] - ds[1L]
par <- c(a1, ak, dk)
names(par) <- c('a1', paste0('ak', 0:(ncat-1)), paste0('d', 0:(ncat-1)))
} else {
stop('traditional2mirt item class not supported', call.=FALSE)
}
par
}
#' Extract parameter variance covariance matrix
#'
#' Extract parameter variance covariance matrix
#'
#' @param object an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{MixedClass}
#'
#' @name vcov-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @export
#' @aliases vcov,SingleGroupClass-method vcov,MixtureClass-method
#' vcov,MultipleGroupClass-method vcov,MixedClass-method vcov,DiscreteClass-method
#' @rdname vcov-method
#' @examples
#'
#' \dontrun{
#' x <- mirt(Science, 1, SE=TRUE)
#' vcov(x)
#'
#' }
setMethod(
f = "vcov",
signature = signature(object = 'SingleGroupClass'),
definition = function(object){
extract.mirt(object, 'vcov')
}
)
#' Extract log-likelihood
#'
#' Extract the observed-data log-likelihood.
#'
#' @param object an object of class \code{SingleGroupClass},
#' \code{MultipleGroupClass}, or \code{MixedClass}
#'
#' @name logLik-method
#' @aliases logLik,SingleGroupClass-method logLik,MixtureClass-method
#' logLik,MultipleGroupClass-method logLik,MixedClass-method logLik,DiscreteClass-method
#' @references
#' Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory
#' Package for the R Environment. \emph{Journal of Statistical Software, 48}(6), 1-29.
#' \doi{10.18637/jss.v048.i06}
#' @rdname logLik-method
#' @examples
#'
#' \dontrun{
#' x <- mirt(Science, 1)
#' logLik(x)
#'
#' }
setMethod(
f = "logLik",
signature = signature(object = 'SingleGroupClass'),
definition = function(object){
extract.mirt(object, 'logLik')
}
)
Try the mirt package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.