Nothing
R/03-estimation.R
In mirt: Multidimensional Item Response Theory
Defines functions
ESTIMATION
ESTIMATION <- function(data, model, group, itemtype = NULL, guess = 0, upper = 1,
invariance = '', pars = NULL, constrain = NULL, key = NULL,
parprior = NULL, mixed.design = NULL, customItems = NULL,
customItemsData = NULL, customGroup = NULL,
GenRandomPars = FALSE, large = FALSE,
survey.weights = NULL, latent.regression = NULL,
gpcm_mats=list(), spline_args=list(), monopoly.k=1,
control = list(), ...)
{
start.time <- proc.time()[3L]
dots <- list(...)
if(!is.null(itemtype))
itemtype <- ifelse(itemtype == 'grsm', 'grsmIRT', itemtype)
if(missing(data)) missingMsg('data')
if(length(unique(colnames(data))) != ncol(data))
stop('items must have unique names in data input', call.=FALSE)
if(missing(model)) missingMsg('model')
if(missing(group)) missingMsg('group')
if(!(is.factor(group) || is.character(group)) || length(group) != nrow(data))
stop('group input provided is not valid', call.=FALSE)
if(is.character(large) && large == 'return'){
Data <- opts <- list()
opts$zeroExtreme <- FALSE
opts$dentype <- 'default'
if(!is.null(dots$technical$zeroExtreme)) opts$zeroExtreme <- dots$technical$zeroExtreme
data <- as.matrix(data)
if(is.numeric(data))
data <- matrix(as.integer(data), nrow(data), dimnames=list(rownames(data), colnames(data)))
rownames(data) <- 1L:nrow(data)
if(is.null(colnames(data)))
colnames(data) <- paste0('Item.', 1L:ncol(data))
Data$data <- data
Data$group <- factor(group)
Data$groupNames <- levels(Data$group)
Data$ngroups <- length(Data$groupNames)
Data$nitems <- ncol(data)
K <- apply(Data$data, 2L, function(x) length(unique(na.omit(x))))
tmp <- list(...)
if(!is.null(tmp$technical$customK)) K <- tmp$technical$customK
itemloc <- c(1L, cumsum(K) + 1L)
Names <- NULL
for(i in 1L:length(K))
Names <- c(Names, paste("Item.",i,"_",1L:K[i],sep=""))
PrepListFull <- list(K=K, itemloc=itemloc,
Names=Names, itemnames=colnames(Data$data))
} else {
if(missing(data) || is.null(nrow(data)))
stop('data argument is required', call.=FALSE)
if(missing(model) || !is(model, 'numeric') && !is(model, 'mirt.model') &&
!is(model, 'character'))
stop('model argument (numeric, character, or from mirt.model() function) is required',
call.=FALSE)
if(!(is.factor(group) || is.character(group)) || length(group) != nrow(data))
stop('group input provided is not valid', call.=FALSE)
if(is.matrix(itemtype)){
itemtypefull <- itemtype
itemtype <- itemtypefull[1L,]
} else itemtypefull <- NULL
if(!is.null(itemtype))
stopifnot(is(itemtype, 'character'))
if(!is.null(constrain))
stopifnot(is(constrain, 'list'))
if(!is.null(parprior))
stopifnot(is(parprior, 'list'))
if(!is.null(customItems))
stopifnot(is(customItems, 'list'))
if(!is.null(customItemsData))
stopifnot(is(customItemsData, 'list') || length(customItemsData) == ncol(data))
if(!is.null(customGroup))
stopifnot(is(customGroup, 'GroupPars') || is.list(customGroup))
stopifnot(is(invariance, 'character'))
stopifnot(is(GenRandomPars, 'logical'))
stopifnot(is(large, 'logical') || is(large, 'list') || is(large, 'character'))
opts <- makeopts(GenRandomPars=GenRandomPars,
hasCustomGroup=!is.null(customGroup), ...)
if(opts$Moptim == 'NR'){
if(is.null(control$tol)) control$tol <- opts$TOL/1000
if(is.null(control$maxit)) control$maxit <- 50L
}
if(opts$Moptim == 'nlminb')
if(is.null(control$rel.tol)) control$rel.tol <- opts$TOL/100
if(opts$dentype == "discrete" && is.null(customGroup)){
den <- Theta_discrete_den
par <- if(!is.null(dots$structure)){
if(is.null(opts$technical$customTheta))
stop('customTheta must be defined when using structure input', call.=FALSE)
opts$structure <- model.matrix(dots$structure, as.data.frame(opts$technical$customTheta))
tmpnfact <- ncol(opts$technical$customTheta)
numeric(ncol(opts$structure))
} else if(is.null(opts$technical$customTheta)){
tmpnfact <- model
numeric(model-1L)
} else {
tmpnfact <- ncol(opts$technical$customTheta)
numeric(nrow(opts$technical$customTheta) - 1L)
}
if(length(par)){
names(par) <- paste0('c', 1:length(par))
est <- rep(TRUE, length(par))
} else {
par <- c(c = 0)
est <- FALSE
}
customGroup <- createGroup(par=par, est=est, den=den, nfact=tmpnfact,
gen=function(object) rnorm(length(object@par), 0, 1/2))
customGroup@itemclass <- -1L
rm(tmpnfact)
}
if(any(itemtype == 'ULL')){
opts$dentype <- 'custom'
den <- function(obj, Theta){
par <- obj@par
dlnorm(Theta, meanlog = par[1L], sdlog = par[2L])
}
opts$theta_lim <- c(.01, opts$theta_lim[2L]^2)
par <- c(meanlog=0, sdlog=1)
est <- c(FALSE, FALSE)
customGroup <- createGroup(par=par, est=est, den=den, nfact=1L,
gen=function(object) rnorm(length(object@par), 0, 1/2))
}
if(!is.null(survey.weights)){
stopifnot(opts$method %in% c('EM', 'QMCEM', 'MCEM'))
stopifnot(length(survey.weights) == nrow(data))
}
if(any(is.na(group))){
if(opts$message)
message('NA values in group removed, along with associated rows in data')
data <- data[!is.na(group), , drop=FALSE]
group <- group[!is.na(group)]
}
if(!is.null(customItems) || any(itemtype %in% Experimental_itemtypes()))
opts$calcNull <- FALSE
opts$times <- list(start.time=start.time)
# on exit, reset the seed to override internal
if(opts$method == 'MHRM' || opts$method == 'MIXED' || opts$method == 'SEM' &&
opts$plausible.draws == 0L)
on.exit(set.seed((as.numeric(Sys.time()) - floor(as.numeric(Sys.time()))) * 1e8))
#change itemtypes if NULL.MODEL
if(opts$NULL.MODEL){
constrain <- NULL
if(!is.null(itemtype)){
itemtype[itemtype == 'grsm'] <- 'graded'
itemtype[itemtype == 'rsm'] <- 'gpcm'
itemtype[itemtype == '3PL' | itemtype == '3PLu' | itemtype == '4PL'] <- '2PL'
itemtype[itemtype == '3PLNRM' | itemtype == '3PLuNRM' | itemtype == '4PLNRM'] <- '2PLNRM'
itemtype[itemtype == 'spline'] <- '2PL'
}
}
if(!is.null(itemtype)){
if(any(itemtype == 'spline') && !(opts$method %in% c('EM', 'QMCEM', 'MCEM')))
stop('spline itemtype only supported for EM algorithm', call.=FALSE)
}
if(length(group) != nrow(data))
stop('length of group not equal to number of rows in data.', call.=FALSE)
if(any(is.na(group)))
stop('Unknown group memberships cannot be estimated. Please remove the NA values in group
and the associated rows in the data input.', call.=FALSE)
#####
### uncomment when building html examples
# opts$verbose <- FALSE
####
opts$times$start.time.Data <- proc.time()[3L]
Data <- list()
data <- as.matrix(data)
if(!all(apply(data, 2L, class) %in% c('integer', 'numeric')))
stop('Input data must be integer or numeric values only', call.=FALSE)
if(is.numeric(data))
data <- matrix(as.integer(data), nrow(data),
dimnames=list(rownames(data), colnames(data)))
rownames(data) <- 1L:nrow(data)
if(is.null(colnames(data)))
colnames(data) <- paste0('Item.', 1L:ncol(data))
if(nrow(data) > 1L && is.null(opts$technical$customK)){
if(!is.null(key)){
key <- sapply(1L:length(key), function(ind, data, key){
if(is.na(key[ind])) return(NA)
s <- sort(unique(data[,ind]))
se <- min(s, na.rm = TRUE):(length(s) + min(s) - 1L)
ret <- se[s == key[ind]]
if(!length(ret)) ret <- NA
ret
}, data=data, key=key)
}
data <- remap.distance(data, message = opts$message)
}
Data$rowID <- 1L:nrow(data)
Data$completely_missing <- which(rowSums(is.na(data)) == ncol(data))
if(length(Data$completely_missing)){
if(opts$warn)
warning('data contains response patterns with only NAs', call.=FALSE)
pick <- rowSums(is.na(data)) != ncol(data)
Data$rowID <- unname(which(pick))
data <- subset(data, pick)
group <- subset(group, pick)
if(!is.null(latent.regression) || !is.null(mixed.design))
covdata <- covdata[-Data$completely_missing, , drop=FALSE]
}
Data$data <- data
if(is.null(opts$grsm.block)) Data$grsm.block <- rep(1L, ncol(data))
else Data$grsm.block <- opts$grsm.block
if(is.null(opts$rsm.block)) Data$rsm.block <- rep(1L, ncol(data))
else Data$rsm.block <- opts$rsm.block
Data$group <- factor(group)
Data$groupNames <- levels(Data$group)
if(any(grepl('-', Data$groupNames)))
stop('Group names cannot contain a dash (-) character', call.=FALSE)
Data$ngroups <- if(!is.null(opts$ngroups)) opts$ngroups else length(Data$groupNames)
if(!is.null(itemtypefull))
if(Data$ngroups != nrow(itemtypefull))
stop("number of groups does not match number of rows in itemtype", call.=FALSE)
Data$nitems <- ncol(Data$data)
Data$N <- nrow(Data$data)
Data$mins <- suppressWarnings(apply(data, 2L, min, na.rm=TRUE))
Data$mins[!is.finite(Data$mins)] <- 0L
if(!is.null(opts$technical$mins)) Data$mins <- opts$technical$mins
if(is.character(model)){
tmp <- any(sapply(colnames(data), grepl, x=model))
model <- mirt.model(model, itemnames = if(tmp) colnames(data) else NULL)
}
oldmodel <- model
PrepList <- vector('list', Data$ngroups)
tmp <- 1L:Data$ngroups
if(opts$dentype == "mixture")
Data$groupNames <- paste0("MIXTURE_", seq_len(Data$ngroups))
names(PrepList) <- Data$groupNames
model <- buildModelSyntax(model, J=Data$nitems, groupNames=Data$groupNames,
itemtype=itemtype)
Data$model <- model
if(!is.null(customGroup)){
if(Data$ngroups == 1L) customGroup <- list(customGroup)
else {
if(Data$ngroups != length(customGroup) && length(customGroup) == 1L){
tmplst <- vector('list', Data$ngroups)
for(g in seq_len(Data$ngroups))
tmplst[[g]] <- customGroup
customGroup <- tmplst
rm(tmplst)
}
}
}
if(!is.null(dots$PrepList)) {
PrepListFull <- PrepList[[1L]] <- dots$PrepList
} else {
PrepListFull <- PrepList[[1L]] <-
PrepData(data=Data$data, model=Data$model, itemtype=itemtype, guess=guess,
upper=upper, parprior=parprior, verbose=opts$verbose,
technical=opts$technical, parnumber=1L, BFACTOR=opts$dentype == 'bfactor',
grsm.block=Data$grsm.block, rsm.block=Data$rsm.block,
mixed.design=mixed.design, customItems=customItems,
customItemsData=customItemsData,
customGroup=customGroup[[1L]], spline_args=spline_args, monopoly.k=monopoly.k,
fulldata=opts$PrepList[[1L]]$fulldata, key=key, opts=opts,
gpcm_mats=gpcm_mats, internal_constraints=opts$internal_constraints,
dcIRT_nphi=opts$dcIRT_nphi, dentype=opts$dentype, item.Q=opts$item.Q)
if(!is.null(dots$Return_PrepList)) return(PrepListFull)
if(!is.null(itemtypefull)){
for(g in 2L:nrow(itemtypefull)){
PrepList[[g]] <-
PrepData(data=Data$data, model=Data$model, itemtype=itemtypefull[g,], guess=guess,
upper=upper, parprior=parprior, verbose=opts$verbose,
technical=opts$technical, parnumber=1L, BFACTOR=opts$dentype == 'bfactor',
grsm.block=Data$grsm.block, rsm.block=Data$rsm.block,
mixed.design=mixed.design, customItems=customItems,
customItemsData=customItemsData,
customGroup=customGroup[[1L]], spline_args=spline_args, monopoly.k=monopoly.k,
fulldata=opts$PrepList[[1L]]$fulldata, key=key, opts=opts,
gpcm_mats=gpcm_mats, internal_constraints=opts$internal_constraints,
dcIRT_nphi=opts$dcIRT_nphi, dentype=opts$dentype, item.Q=opts$item.Q)
}
}
}
if(!is.null(opts$structure)){
PrepList[[1L]]$pars[[Data$nitems+1L]]@structure <- opts$structure
PrepList[[1L]]$pars[[Data$nitems+1L]]@parnames <-
names(PrepList[[1L]]$pars[[Data$nitems+1L]]@est) <- colnames(opts$structure)
PrepList[[1L]]$pars[[Data$nitems+1L]]@itemclass <- -999L
}
parnumber <- max(PrepList[[1L]]$pars[[Data$nitems+1L]]@parnum) + 1L
attr(PrepListFull$pars, 'nclasspars') <- attr(PrepList[[1L]]$pars, 'nclasspars') <-
matrix(sapply(PrepListFull$pars, function(y) length(y@parnum)), nrow=1L)
for(g in seq_len(Data$ngroups)){
if(g != 1L){
if(is.null(itemtypefull))
PrepList[[g]] <- list(pars=PrepList[[1L]]$pars)
else attr(PrepList[[g]]$pars, 'nclasspars') <-
sapply(PrepList[[g]]$pars, function(y) length(y@parnum))
for(i in 1L:length(PrepList[[g]]$pars)){
PrepList[[g]]$pars[[i]]@parnum <- parnumber:(parnumber +
length(PrepList[[g]]$pars[[i]]@parnum) - 1L)
parnumber <- max(PrepList[[g]]$pars[[i]]@parnum) + 1L
}
if(!is.null(customGroup)){
customGroup[[g]]@parnum <- PrepList[[g]]$pars[[Data$nitems + 1L]]@parnum
PrepList[[g]]$pars[[Data$nitems + 1L]] <- customGroup[[g]]
}
}
}
if(length(mixed.design$random) > 0L){
for(i in seq_len(length(mixed.design$random))){
mixed.design$random[[i]]@parnum <- parnumber:(parnumber - 1L +
length(mixed.design$random[[i]]@par))
parnumber <- max(mixed.design$random[[i]]@parnum) + 1L
}
}
if(!is.null(latent.regression)){
if(length(PrepListFull$prodlist))
stop('Polynomial combinations currently not supported when latent regression effects are used', call.=FALSE)
lrPars <- make.lrdesign(df=latent.regression$df, formula=latent.regression$formula,
factorNames=PrepListFull$factorNames, EM=latent.regression$EM,
TOL=opts$TOL)
lrPars@parnum <- parnumber:(parnumber - 1L + length(lrPars@par))
parnumber <- max(lrPars@parnum) + 1L
if(opts$dentype == 'discrete'){
tmp <- matrix(1L:length(lrPars@beta), nrow(lrPars@beta), ncol(lrPars@beta))
tmp2 <- tmp[1, ]
lrPars@est[tmp2[-length(tmp2)]] <- TRUE
lrPars@est[tmp[,ncol(tmp)]] <- FALSE
PrepList$all$pars[[ncol(data) + 1L]]@est <- FALSE
}
} else lrPars <- list()
if(length(latent.regression$lr.random) > 0L){
for(i in seq_len(length(latent.regression$lr.random))){
latent.regression$lr.random[[i]]@parnum <- parnumber:(parnumber - 1L +
length(latent.regression$lr.random[[i]]@par))
parnumber <- max(latent.regression$lr.random[[i]]@parnum) + 1L
}
}
}
if(!is.null(opts$PrepList)){
PrepList <- opts$PrepList
} else {
if(!is.list(large)){
tmpdata <- Data$data
if(opts$zeroExtreme){
n_is_na <- rowSums(is.na(tmpdata))
sums <- colSums(t(tmpdata) - Data$mins, na.rm=TRUE)
maxs <- apply(tmpdata, 2L, function(x) length(unique(na.omit(x)))) - 1L
if(is.null(survey.weights)) survey.weights <- rep(1, nrow(tmpdata))
survey.weights[sums == 0L | (sums + n_is_na) == sum(maxs)] <- 0
}
tmptabdata <- if(is.logical(large) && large){
maketabDataLarge(tmpdata=tmpdata, group=Data$group,
groupNames=if(opts$dentype != 'mixture') Data$groupNames else 'full',
nitem=Data$nitems, K=PrepListFull$K, itemloc=PrepListFull$itemloc,
Names=PrepListFull$Names, itemnames=PrepListFull$itemnames,
survey.weights=survey.weights)
} else {
maketabData(tmpdata=tmpdata, group=Data$group,
groupNames=if(opts$dentype != 'mixture') Data$groupNames else 'full',
nitem=Data$nitems, K=PrepListFull$K, itemloc=PrepListFull$itemloc,
Names=PrepListFull$Names, itemnames=PrepListFull$itemnames,
survey.weights=survey.weights)
}
if(is.character(large) && large == 'return'){
return(tmptabdata)
} else large <- tmptabdata
rm(tmpdata, tmptabdata)
}
Data$tabdatalong <- large$tabdata
Data$tabdata <- large$tabdata2
for(g in seq_len(Data$ngroups)){
if(opts$dentype == 'mixture'){
Data$fulldata[[g]] <- PrepListFull$fulldata
Data$Freq[[g]] <- large$Freq[[1L]]
} else {
Data$fulldata[[g]] <- PrepListFull$fulldata[Data$group == Data$groupNames[g],
, drop=FALSE]
Data$Freq[[g]] <- large$Freq[[g]]
}
}
}
if(opts$returnPrepList) return(PrepList)
if(opts$dentype == 'bfactor'){
#better start values
if((PrepList[[1L]]$nfact - attr(model, 'nspec')) == 1L){
nfact <- PrepListFull$nfact
for(g in seq_len(Data$ngroups)){
for(i in seq_len(Data$nitems)){
tmp <- PrepList[[g]]$pars[[i]]@par[1L:nfact]
tmp2 <- tmp[1L]
tmp[PrepList[[g]]$pars[[i]]@est[1L:nfact]] <-
tmp[PrepList[[g]]$pars[[i]]@est[1L:nfact]]/2
tmp[1L] <- tmp2
PrepList[[g]]$pars[[i]]@par[1L:nfact] <- tmp
}
}
}
}
PrepList <- UpdateParameters(PrepList, model, groupNames=Data$groupNames)
if(GenRandomPars){
for(g in seq_len(Data$ngroups))
for(i in seq_len(length(PrepList[[g]]$pars)))
PrepList[[g]]$pars[[i]] <- GenRandomPars(PrepList[[g]]$pars[[i]])
}
if(opts$dentype == "discrete"){
PrepList[[1L]]$exploratory <- FALSE
if(is.null(opts$technical$customTheta))
opts$technical$customTheta <- diag(PrepList[[1L]]$nfact)
}
RETURNVALUES <- FALSE
SUPPLIED_STARTS <- FALSE
if(!is.null(pars)){
if(is(pars, 'data.frame')){
SUPPLIED_STARTS <- TRUE
PrepList <- UpdatePrepList(PrepList, pars, random=mixed.design$random,
lrPars=lrPars, lr.random=latent.regression$lr.random,
MG = TRUE)
mixed.design$random <- attr(PrepList, 'random')
latent.regression$lr.random <- attr(PrepList, 'lr.random')
if(any(pars$class == 'lrPars')) lrPars <- update.lrPars(pars, lrPars)
attr(PrepList, 'random') <- NULL
attr(PrepList, 'lr.random') <- NULL
}
if(!is.null(attr(pars, 'values')) || (is.character(pars)))
RETURNVALUES <- TRUE
}
pars <- vector('list', Data$ngroups)
for(g in seq_len(Data$ngroups))
pars[[g]] <- PrepList[[g]]$pars
nitems <- Data$nitems
Data$K <- PrepList[[1L]]$K
nfact <- PrepList[[1L]]$pars[[nitems+1L]]@nfact
if(nfact != 1L && any(c('Rasch') %in% itemtype ) && PrepList[[1L]]$exploratory)
stop('Rasch itemtype is for confirmatory models only.', call.=FALSE)
nLambdas <- PrepList[[1L]]$pars[[1L]]@nfact
if(is.null(constrain)) constrain <- list()
nspec <- ifelse(!is.null(attr(model, 'nspec')), attr(model, 'nspec'), 1L)
#default MG uses configural model (independent groups but each identified)
if('free_means' %in% invariance ){ #Free factor means (means 0 for ref)
if(opts$dentype == 'bfactor'){
for(g in 2L:Data$ngroups)
pars[[g]][[nitems + 1L]]@est[1L:(nfact-nspec)] <- TRUE
} else {
for(g in 2L:Data$ngroups)
pars[[g]][[nitems + 1L]]@est[1L:nfact] <- TRUE
}
}
dummymat <- matrix(FALSE, pars[[1L]][[nitems + 1L]]@nfact, pars[[1L]][[nitems + 1L]]@nfact)
if(any('free_means' %in% invariance)){ #Free means
if(all(sapply(PrepList[[1]]$pars, function(x) class(x)) %in%
c(ordinal_itemtypes(), 'GroupPars'))){
TS <- rowMeans(Data$data, na.rm=TRUE)
gmus <- tapply(TS, Data$group, mean, na.rm=TRUE)
gsd <- sd(TS[Data$group == Data$groupNames[1L]], na.rm=TRUE)
gmuscaled <- (gmus - gmus[Data$groupNames[1L]]) / gsd
if(all(is.finite(gmuscaled))){
for(i in 1L:Data$ngroups){
tmp <- pars[[i]][[Data$nitems+1L]]
tmp@par[tmp@est & grepl('MEAN_', names(tmp@est))] <- gmuscaled[i]
pars[[i]][[Data$nitems+1L]] <- tmp
}
}
}
}
if(any(c('free_var', 'free_vars') %in% invariance)){ #Free factor vars (vars 1 for ref)
if(opts$dentype == 'bfactor'){
tmp <- dummymat[1L:(nfact-nspec),1L:(nfact-nspec), drop=FALSE]
diag(tmp) <- TRUE
dummymat[1L:(nfact-nspec),1L:(nfact-nspec)] <- tmp
} else {
diag(dummymat) <- TRUE
}
tmp <- dummymat[lower.tri(dummymat, TRUE)]
if(opts$dentype %in% c('Davidian', 'EHW')){
for(g in 2L:Data$ngroups)
pars[[g]][[nitems + 1L]]@est[2L] <- TRUE
} else {
for(g in 2L:Data$ngroups){
pars[[g]][[nitems + 1L]]@est <- pars[[g]][[nitems + 1L]]@est |
c(pars[[g]][[nitems + 1L]]@est[1L:pars[[g]][[nitems + 1L]]@nfact], tmp)
names(pars[[g]][[nitems + 1L]]@est) <- names(pars[[g]][[nitems + 1L]]@par)
pars[[g]][[nitems + 1L]]@parnames <- names(pars[[g]][[nitems + 1L]]@est)
}
}
}
if(opts$dentype == 'mixture' && !SUPPLIED_STARTS){
tmp <- length(pars[[1L]][[nitems + 1L]]@par)
pars[[1L]][[nitems + 1L]]@est[tmp] <- FALSE
for(g in 1L:Data$ngroups)
pars[[g]][[nitems + 1L]]@par[tmp] <- g - 1
names(PrepList) <- Data$groupNames
}
constrain <- UpdateConstrain(pars=pars, constrain=constrain, invariance=invariance, nfact=Data$nfact,
nLambdas=nLambdas, J=nitems, ngroups=Data$ngroups, PrepList=PrepList,
method=opts$method, itemnames=PrepList[[1L]]$itemnames, model=model,
groupNames=Data$groupNames)
pars <- resetPriorConstrain(pars=pars, constrain=constrain,
nconstrain=opts$technical$nconstrain)
if(RETURNVALUES){
for(g in seq_len(Data$ngroups))
PrepList[[g]]$pars <- pars[[g]]
return(ReturnPars(PrepList, PrepList[[1L]]$itemnames, lr.random=latent.regression$lr.random,
random=mixed.design$random, lrPars=lrPars, MG = TRUE))
}
startlongpars <- c()
if(opts$NULL.MODEL){
constrain <- list()
for(g in seq_len(length(pars))){
for(i in seq_len(nitems)){
pars[[g]][[i]] <- set_null_model(pars[[g]][[i]])
}
}
}
rr <- 0L
for(g in seq_len(Data$ngroups)){
if(g > 1L && opts$dentype == 'mixture') break
r <- Data$Freq[[g]]
rr <- rr + r
}
df <- if(opts$dentype == 'mixture') prod(PrepList[[1L]]$K) - 1
else Data$ngroups * (prod(PrepList[[1L]]$K) - 1)
if(df > 1e10) df <- 1e10
nestpars <- nconstr <- 0L
for(g in seq_len(Data$ngroups))
for(i in seq_len(nitems+1L))
nestpars <- nestpars + sum(pars[[g]][[i]]@est)
if(!is.null(mixed.design$random)){
for(i in seq_len(length(mixed.design$random)))
nestpars <- nestpars + sum(mixed.design$random[[i]]@est)
}
if(length(lrPars))
nestpars <- nestpars + sum(lrPars@est)
if(!is.null(dots$structure)) nestpars <- nestpars - 1L
if(!is.null(opts$technical$nconstrain))
constrain <- c(constrain, opts$technical$nconstrain)
if(length(constrain) > 0L)
for(i in seq_len(length(constrain)))
nconstr <- nconstr + length(constrain[[i]]) - 1L
if(Data$ngroups > 1L && !length(constrain)){
if(opts$warn && any(invariance %in% c('free_means', 'free_var')))
warning('Multiple-group model may not be identified without providing anchor items',
call.=FALSE)
for(j in seq_len(Data$ngroups)){
if(opts$dentype != 'mixture')
tmp <- apply(subset(Data$data, Data$group == Data$groupNames[j]), 2L,
function(x) length(unique(na.omit(x)))) == Data$K
for(i in which(!tmp)){
if(any(PrepList[[j]]$pars[[i]]@est))
stop(paste0('Multiple Group model will not be identified without ',
'proper constraints (groups contain missing data patterns ',
'where item responses have been completely omitted or, alternatively, ',
'the number of categories within each group is not equal to the ',
'total number of categories)'),
call. = FALSE)
}
}
}
nmissingtabdata <- sum(is.na(rowSums(Data$tabdata)))
dfsubtr <- nestpars - nconstr
if(opts$dentype == 'EH') dfsubtr <- dfsubtr + (opts$quadpts - 1L) * Data$ngroups
else if(opts$dentype == 'EHW') dfsubtr <- dfsubtr + (opts$quadpts - 3L) * Data$ngroups
if(df <= dfsubtr)
stop('Too few degrees of freedom. There are only ', df, ' degrees of freedom but ',
dfsubtr, ' parameters were freely estimated.', call.=FALSE)
df <- df - dfsubtr
if(!is.null(customItems)){
for(g in seq_len(Data$ngroups))
PrepList[[g]]$exploratory <- FALSE
}
G2group <- numeric(Data$ngroups)
DERIV <- vector('list', Data$ngroups)
for(g in seq_len(Data$ngroups)){
DERIV[[g]] <- vector('list', Data$nitems)
for(i in seq_len(Data$nitems))
DERIV[[g]][[i]] <- selectMethod(Deriv, c(class(pars[[g]][[i]]), 'matrix'))
}
Ls <- makeLmats(pars, constrain, random = mixed.design$random,
lr.random=latent.regression$lr.random, lrPars=lrPars,
nconstrain=opts$technical$nconstrain)
CUSTOM.IND <- which(sapply(pars[[1L]], class) %in% Use_R_ProbTrace())
SLOW.IND <- which(sapply(pars[[1L]], class) %in% Use_R_Deriv())
if(pars[[1]][[length(pars[[1L]])]]@itemclass %in% c(-1L, -999L))
SLOW.IND <- c(SLOW.IND, length(pars[[1L]]))
if(opts$dentype != 'Gaussian' && opts$method %in% c('MHRM', 'MIXED', 'SEM'))
stop('Non-Gaussian densities not currently supported with MHRM algorithm')
#warnings
wmsg <- 'Lower and upper bound parameters (g and u) should use \'norm\' (i.e., logit) prior'
for(g in seq_len(length(pars))){
for(i in seq_len(length(pars[[1L]]))){
if(is(pars[[g]][[i]], 'dich')){
pt <- pars[[g]][[i]]@prior.type
if(!(pt[length(pt)-1L] %in% c(0L, 1L, 4L))) warning(wmsg, call.=FALSE)
if(!(pt[length(pt)] %in% c(0L, 1L, 4L))) warning(wmsg, call.=FALSE)
next
} else if(is(pars[[g]][[i]], 'partcomp')){
pt <- pars[[g]][[i]]@prior.type
if(!(pt[length(pt)-1L] %in% c(0L, 1L))) warning(wmsg, call.=FALSE)
if(!(pt[length(pt)] %in% c(0L, 1L))) warning(wmsg, call.=FALSE)
next
} else if(is(pars[[g]][[i]], 'nestlogit')){
pt <- pars[[g]][[i]]@prior.type
if(!(pt[nfact + 2L] %in% c(0L, 1L))) warning(wmsg, call.=FALSE)
if(!(pt[nfact + 3L] %in% c(0L, 1L))) warning(wmsg, call.=FALSE)
next
}
}
}
SEMconv <- NA
Data$wmiss <- if(length(lrPars)) with(Data, 1/rowMeans(!is.na(data)) / nitems)
else with(Data, 1/rowMeans(!is.na(tabdata)) / nitems)
opts$times$end.time.Data <- proc.time()[3L]
#EM estimation
opts$times$start.time.Estimate <- proc.time()[3L]
logLik <- G2 <- SElogLik <- NaN
logPrior <- 0
Pl <- vector('list', Data$ngroups)
if(opts$method %in% c('EM', 'BL', 'QMCEM', 'MCEM')){
logLik <- G2 <- SElogLik <- 0
if(length(lrPars)){
if(opts$SE && !(opts$SE.type %in% c('complete', 'forward', 'central', 'Richardson')))
stop('Information matrix method for latent regression estimates not supported',
call.=FALSE)
opts$full <- TRUE
} else opts$full <- FALSE
temp <- matrix(0L,nrow=nitems,ncol=nspec)
sitems <- matrix(0L, nrow=sum(PrepList[[1L]]$K), ncol=nspec)
specific <- NULL
if(opts$dentype == "discrete"){
theta <- 0
Theta <- opts$technical$customTheta
opts$quadpts <- nrow(Theta)
if(pars[[1L]][[1L]]@nfact != ncol(Theta))
stop("mirt.model definition does not have same number of traits/attributes as customTheta input", call.=FALSE)
} else {
if(is.null(opts$quadpts)){
tmp <- if(opts$dentype == 'bfactor') PrepList[[1L]]$nfact - attr(model, 'nspec') + 1L
else nfact
opts$quadpts <- select_quadpts(tmp)
}
if(opts$quadpts < 3 && opts$warn) warning('Should use more than 2 quadpts', call.=FALSE)
theta <- 1
if(!(opts$method %in% c('QMCEM', 'MCEM')))
theta <- as.matrix(seq(opts$theta_lim[1L], opts$theta_lim[2L],
length.out = opts$quadpts))
if(opts$dentype == 'bfactor'){
specific <- attr(oldmodel, 'specific')
specific[is.na(specific)] <- 1L
for(i in seq_len(nitems)) temp[i, specific[i]] <- 1L
ind <- 1L
for(i in seq_len(nitems)){
for(j in seq_len(PrepList[[1L]]$K[i])){
sitems[ind, ] <- temp[i, ]
ind <- ind + 1L
}
}
nfact2 <- PrepList[[1L]]$nfact - attr(model, 'nspec') + 1L
Theta <- thetaComb(theta, nfact2)
Theta <- cbind(Theta[,1L:(nfact2-1L),drop=FALSE],
matrix(Theta[,nfact2], nrow=nrow(Theta), ncol=ncol(sitems)))
} else {
if(opts$method %in% c('QMCEM', 'MCEM')){
Theta <- NULL
} else {
if(opts$quadpts^nfact <= opts$MAXQUAD){
if(is.null(opts$technical$customTheta))
Theta <- thetaComb(theta, nfact)
} else stop('Greater than ', opts$MAXQUAD, ' quadrature points.', call.=FALSE)
if(opts$message && nfact > 3L && !(opts$odentype %in% c('custom', 'discrete')))
message('EM quadrature for high dimensional models are better handled
\twith the \"QMCEM\" or \"MCEM\" method')
}
}
if(!is.null(opts$technical$customTheta)){
Theta <- opts$technical$customTheta
if(!is.matrix(Theta)) stop('customTheta input must be a matrix', call.=FALSE)
opts$quadpts <- nrow(Theta)
}
pars <- loadSplinePars(pars, Theta)
} #end Theta def
ESTIMATE <- EM.group(pars=pars, constrain=constrain, Ls=Ls, PrepList=PrepList, Data=Data,
list = list(NCYCLES=opts$NCYCLES, TOL=opts$TOL, MSTEPTOL=opts$MSTEPTOL,
nfactNames=PrepList[[1L]]$nfactNames, theta=theta,
itemloc=PrepList[[1L]]$itemloc, dentype=opts$dentype,
sitems=sitems, specific=specific, NULL.MODEL=opts$NULL.MODEL,
nfact=nfact, constrain=constrain, verbose=opts$verbose,
SE = opts$SE, SE.type=opts$SE.type, delta=opts$delta, quadpts=opts$quadpts,
accelerate=opts$accelerate, CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND,
customPriorFun=opts$customPriorFun, Moptim=opts$Moptim, warn=opts$warn,
message=opts$message, method=opts$method, full=opts$full,
lrPars=lrPars, SE=opts$SE && opts$SE.type == 'numerical', Etable=opts$Etable,
NULL.MODEL=opts$NULL.MODEL, PLCI=opts$PLCI, Norder=opts$Norder,
keep_vcov_PD=opts$keep_vcov_PD, symmetric=opts$technical$symmetric,
MCEM_draws=opts$MCEM_draws, omp_threads=opts$omp_threads),
Theta=Theta, DERIV=DERIV, solnp_args=opts$solnp_args, control=control,
nconstrain=opts$technical$nconstrain)
if(opts$method == 'MCEM')
opts$quadpts <- opts$MCEM_draws(ESTIMATE$cycles)
opts$Moptim <- ESTIMATE$Moptim
lrPars <- ESTIMATE$lrPars
startlongpars <- ESTIMATE$longpars
rlist <- ESTIMATE$rlist
logPrior <- ESTIMATE$logPrior
for(g in seq_len(Data$ngroups)){
if(g > 1L && opts$dentype == 'mixture') break
Pl[[g]] <- rlist[[g]]$expected
Pltmp <- Pl[[g]]
if(opts$full){
rg <- 1
G2group[g] <- NaN
} else {
rg <- Data$Freq[[g]]
Pltmp <- Pltmp[rg != 0]
rg <- rg[rg != 0]
Ng <- sum(rg)
G2group[g] <- 2 * sum(rg * log(rg/(Ng*Pltmp)))
}
G2 <- G2 + G2group[g]
logLik <- logLik + sum(rg*log(Pltmp))
}
} else if(opts$method %in% c('MHRM', 'SEM')){ #MHRM estimation
Theta <- matrix(0, Data$N, nitems)
if(opts$method == 'SEM') opts$NCYCLES <- NA
ESTIMATE <- MHRM.group(pars=pars, constrain=constrain, Ls=Ls, PrepList=PrepList, Data=Data,
list = list(NCYCLES=opts$NCYCLES, BURNIN=opts$BURNIN,
SEMCYCLES=opts$SEMCYCLES, gain=opts$gain,
KDRAWS=opts$KDRAWS, MHDRAWS=opts$MHDRAWS,
TOL=opts$TOL, SE=FALSE, SE.type = 'none',
nfactNames=PrepList[[1L]]$nfactNames,
itemloc=PrepList[[1L]]$itemloc,
nfact=nfact, constrain=constrain, verbose=opts$verbose,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND,
startlongpars=startlongpars,
cand.t.var=opts$technical$MHcand, warn=opts$warn,
message=opts$message, expl=PrepList[[1L]]$exploratory,
plausible.draws=opts$plausible.draws,
MSTEPTOL=opts$MSTEPTOL, Moptim=opts$Moptim,
keep_vcov_PD=opts$keep_vcov_PD),
DERIV=DERIV, solnp_args=opts$solnp_args, control=control)
if(opts$plausible.draws != 0) return(ESTIMATE)
if(opts$SE && (ESTIMATE$converge || !opts$info_if_converged)){
if(opts$verbose)
cat('\nCalculating information matrix...\n')
tmp <- MHRM.group(pars=ESTIMATE$pars, constrain=constrain, Ls=Ls, PrepList=PrepList, Data=Data,
list = list(NCYCLES=opts$MHRM_SE_draws, BURNIN=1L,
SEMCYCLES=opts$SEMCYCLES, gain=opts$gain,
KDRAWS=opts$KDRAWS, MHDRAWS=opts$MHDRAWS,
TOL=opts$SEtol, SE=TRUE, SE.type=opts$SE.type,
nfactNames=PrepList[[1L]]$nfactNames,
itemloc=PrepList[[1L]]$itemloc,
nfact=nfact, constrain=constrain, verbose=FALSE,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND,
startlongpars=ESTIMATE$longpars, plausible.draws=0L,
cand.t.var=opts$technical$MHcand, warn=opts$warn,
message=opts$message, expl=PrepList[[1L]]$exploratory,
MSTEPTOL=opts$MSTEPTOL, Moptim='NR1',
keep_vcov_PD=opts$keep_vcov_PD),
DERIV=DERIV, solnp_args=opts$solnp_args, control=control)
ESTIMATE$pars <- tmp$pars
ESTIMATE$info <- tmp$info
ESTIMATE$fail_invert_info <- tmp$fail_invert_info
ESTIMATE$time <- c(ESTIMATE$time, SE=sum(tmp$time))
rm(tmp)
}
rlist <- vector('list', Data$ngroups)
for(g in seq_len(Data$ngroups))
rlist[[g]]$expected = numeric(1L)
} else if(opts$method == 'MIXED'){
if(is.null(opts$technical$RANDSTART)) opts$technical$RANDSTART <- 100L
if(is.null(opts$technical$BURNIN) && length(mixed.design$random)) opts$BURNIN <- 200L
Theta <- matrix(0, Data$N, nitems)
ESTIMATE <- MHRM.group(pars=pars, constrain=constrain, Ls=Ls,
PrepList=PrepList, random=mixed.design$random, Data=Data,
lrPars=lrPars, lr.random=latent.regression$lr.random,
list = list(NCYCLES=opts$NCYCLES, BURNIN=opts$BURNIN,
SEMCYCLES=opts$SEMCYCLES, gain=opts$gain,
KDRAWS=opts$KDRAWS, MHDRAWS=opts$MHDRAWS,
TOL=opts$TOL, SE.type = 'none',
nfactNames=PrepList[[1L]]$nfactNames,
itemloc=PrepList[[1L]]$itemloc,
nfact=nfact, constrain=constrain, verbose=opts$verbose,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND,
startlongpars=startlongpars, SE=FALSE,
cand.t.var=opts$technical$MHcand, warn=opts$warn,
message=opts$message, expl=FALSE, plausible.draws=0L,
RANDSTART=opts$technical$RANDSTART,
MSTEPTOL=opts$MSTEPTOL, Moptim=opts$Moptim,
keep_vcov_PD=opts$keep_vcov_PD),
DERIV=DERIV, solnp_args=opts$solnp_args, control=control)
if(opts$SE && (ESTIMATE$converge || !opts$info_if_converged)){
if(opts$verbose)
cat('\nCalculating information matrix...\n')
tmp <- MHRM.group(pars=ESTIMATE$pars, constrain=constrain, Ls=Ls,
PrepList=PrepList, random=mixed.design$random, Data=Data,
lrPars=ESTIMATE$lrPars, lr.random=latent.regression$lr.random,
list = list(NCYCLES=opts$MHRM_SE_draws, BURNIN=1L,
SEMCYCLES=opts$SEMCYCLES, gain=opts$gain,
KDRAWS=opts$KDRAWS, MHDRAWS=opts$MHDRAWS,
TOL=opts$SEtol, SE=TRUE, SE.type=opts$SE.type,
nfactNames=PrepList[[1L]]$nfactNames,
itemloc=PrepList[[1L]]$itemloc,
nfact=nfact, constrain=constrain, verbose=FALSE,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND,
startlongpars=ESTIMATE$longpars, plausible.draws=0L,
cand.t.var=opts$technical$MHcand, warn=opts$warn,
message=opts$message, expl=FALSE,
RANDSTART=1L,
MSTEPTOL=opts$MSTEPTOL, Moptim='NR1',
keep_vcov_PD=opts$keep_vcov_PD),
DERIV=DERIV, solnp_args=opts$solnp_args, control=control)
ESTIMATE$pars <- tmp$pars
ESTIMATE$random <- tmp$random
ESTIMATE$lrPars <- tmp$lrPars
ESTIMATE$lr.random <- tmp$lr.random
ESTIMATE$info <- tmp$info
ESTIMATE$fail_invert_info <- tmp$fail_invert_info
ESTIMATE$time <- c(ESTIMATE$time, SE=sum(tmp$time))
rm(tmp)
}
rlist <- vector('list', Data$ngroups)
for(g in seq_len(Data$ngroups))
rlist[[g]]$expected = numeric(1L)
}
for(g in seq_len(length(pars))){
for(i in seq_len(length(pars[[1L]]))){
if(is(pars[[g]][[i]], 'dich')){
tmp <- ESTIMATE$pars[[g]][[i]]@par
nms <- ESTIMATE$pars[[g]][[i]]@parnames
if(tmp[nms == 'g'] > tmp[nms == 'u']){
if(opts$warn)
warning('g parameter greater than u detected. Model did not converge', call.=FALSE)
ESTIMATE$converge <- FALSE
}
}
}
}
opts$times$end.time.Estimate <- proc.time()[3L]
if(opts$logLik_if_converged && !ESTIMATE$converge) opts$draws <- 0
opts$times$start.time.SE <- proc.time()[3L]
if(!opts$NULL.MODEL && opts$SE){
tmp <- ESTIMATE
if(opts$verbose && !(opts$method %in% c('MHRM', 'MIXED', 'SEM')))
cat('\n\nCalculating information matrix...\n')
if(opts$SE.type %in% c('complete', 'Oakes') && opts$method %in% c('EM', 'QMCEM')){
opts$times$start.time.SE <- ESTIMATE$start.time.SE
ESTIMATE <- loadESTIMATEinfo(info=-ESTIMATE$hess, ESTIMATE=ESTIMATE, constrain=constrain,
warn=opts$warn)
} else if(opts$SE.type == 'SEM' && opts$method == 'EM'){
collectLL <- as.numeric(ESTIMATE$collectLL)
collectLL <- exp(c(NA, collectLL) - c(collectLL, NA))
from <- suppressWarnings(max(which(collectLL <= opts$SEM_from)))
if(from < 1L) from <- 1L
to <- min(which(collectLL >= opts$SEM_to))
dontrun <- FALSE
if(from == to){
if(opts$warn)
warning('SEM window is too small to compute information matrix.
Consider changing the starting values', call.=FALSE)
dontrun <- TRUE
}
lengthsplit <- do.call(c, lapply(strsplit(names(ESTIMATE$correct), 'COV_'), length))
lengthsplit <- lengthsplit + do.call(c, lapply(strsplit(names(ESTIMATE$correct), 'MEAN_'), length))
is.latent <- lengthsplit > 2L
if(!dontrun){
if(ESTIMATE$cycles <= 10L)
if(opts$warn)
warning('Very few EM cycles performed. Consider decreasing TOL further to
increase EM iteration count or starting farther away from ML estimates by
passing the \'GenRandomPars = TRUE\' argument')
estmat <- matrix(FALSE, length(ESTIMATE$correction), length(ESTIMATE$correction))
DM <- estmat + 0
diag(estmat) <- TRUE
if(!opts$technical$parallel){
ncores <- .mirtClusterEnv$ncores
.mirtClusterEnv$ncores <- 1L
}
DM <- myLapply(1L:ncol(estmat), FUN=SE.SEM, estmat=estmat, pars=ESTIMATE$pars, constrain=constrain, Data=Data,
list = list(NCYCLES=opts$NCYCLES, TOL=opts$SEtol, MSTEPTOL=opts$MSTEPTOL,
nfactNames=PrepList[[1L]]$nfactNames, theta=theta,
itemloc=PrepList[[1L]]$itemloc, keep_vcov_PD=opts$keep_vcov_PD,
sitems=sitems, specific=specific, NULL.MODEL=opts$NULL.MODEL,
nfact=nfact, constrain=constrain, verbose=opts$verbose,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND, Moptim=ESTIMATE$Moptim,
EHPrior=ESTIMATE$Prior, warn=opts$warn, dentype=opts$dentype,
message=opts$message, full=opts$full, lrPars=lrPars, method=opts$method),
Theta=Theta, theta=theta, ESTIMATE=ESTIMATE, from=from, to=to,
DERIV=DERIV, is.latent=is.latent, Ls=Ls, PrepList=PrepList,
solnp_args=opts$solnp_args, control=control)
SEMconv <- sapply(DM, function(x) all(attr(x, 'converged')))
if(!all(SEMconv)){
warning(sprintf(c('%i parameters did not converge in numerical SEM derivative.\n',
'Try using different starting values or passing GenRandomPars=TRUE'),
sum(!SEMconv)),
call.=FALSE)
SEMconv <- FALSE
} else SEMconv <- TRUE
DM <- do.call(rbind, DM)
if(!opts$technical$parallel)
.mirtClusterEnv$ncores <- ncores
ESTIMATE$pars <- reloadPars(longpars=ESTIMATE$longpars, pars=ESTIMATE$pars,
ngroups=Data$ngroups, J=Data$nitems)
DM[, is.latent] <- DM[is.latent, ]
DM[is.latent, is.latent] <- 0
info <- try(solve(-solve(ESTIMATE$hess) %*% solve(diag(ncol(DM)) - DM)), silent=TRUE)
info[,is.latent] <- t(info[is.latent, ,drop=FALSE])
if(opts$technical$symmetric) info <- (info + t(info)) / 2
if(is(info, 'try-error')){
if(opts$warn)
warning('Information matrix in SEM could not be computed due to instability',
call.=FALSE)
} else ESTIMATE <- loadESTIMATEinfo(info=info, ESTIMATE=ESTIMATE, constrain=constrain,
warn=opts$warn)
}
} else if(opts$SE.type == 'numerical' && opts$method == 'BL'){
ESTIMATE <- loadESTIMATEinfo(info=-ESTIMATE$hess, ESTIMATE=ESTIMATE, constrain=constrain,
warn=opts$warn)
} else if(opts$SE.type %in% c('Richardson', 'forward', 'central') &&
!(opts$method %in% c('MHRM', 'SEM', 'MIXED'))){
ESTIMATE <- SE.Numerical(pars=ESTIMATE$pars, Theta=ESTIMATE$Theta, theta=theta, PrepList=PrepList, Data=Data,
dentype=opts$dentype, itemloc=PrepList[[1L]]$itemloc, ESTIMATE=ESTIMATE,
constrain=constrain, Ls=Ls, specific=oldmodel, sitems=sitems,
CUSTOM.IND=CUSTOM.IND, EHPrior=ESTIMATE$Prior, warn=opts$warn, type=opts$SE.type,
delta=opts$delta, lrPars=ESTIMATE$lrPars, omp_threads=opts$omp_threads)
} else if(opts$SE.type %in% c('MHRM', 'FMHRM') && opts$method == 'EM'){
if(opts$dentype %in% c('EH', 'EHW'))
stop('MHRM standard error methods not available when using empirical histograms', call.=FALSE)
ESTIMATE <- MHRM.group(pars=pars, constrain=constrain, Ls=Ls, PrepList=PrepList, Data=Data,
list = list(NCYCLES=1000L, BURNIN=1L, SEMCYCLES=opts$SEMCYCLES,
KDRAWS=opts$KDRAWS, MHDRAWS=opts$MHDRAWS,
TOL=opts$SEtol, SE=TRUE, SE.type=opts$SE.type,
gain=opts$gain, nfactNames=PrepList[[1L]]$nfactNames,
itemloc=PrepList[[1L]]$itemloc,
nfact=nfact, constrain=constrain, verbose=FALSE, expl=FALSE,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND, message=opts$message,
startlongpars=startlongpars, SE=opts$SE, warn=opts$warn,
plausible.draws=0L, MSTEPTOL=opts$MSTEPTOL, Moptim='NR1',
keep_vcov_PD=opts$keep_vcov_PD),
DERIV=DERIV, solnp_args=opts$solnp_args, control=control)
} else if(any(opts$SE.type %in% c('crossprod', 'Louis', 'sandwich.Louis', 'sandwich')) &&
!(opts$method %in% c('MHRM', 'SEM', 'MIXED'))){
if(opts$method %in% c('QMCEM', 'MCEM')){
if(opts$warn)
warning(sprintf('SE.type not supported when using method = \"%s\"', opts$method),
call. = FALSE)
} else {
ESTIMATE <- SE.simple(PrepList=PrepList, ESTIMATE=ESTIMATE, Theta=Theta, Data=Data,
constrain=constrain, Ls=Ls, N=nrow(data), type=opts$SE.type,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND, warn=opts$warn,
message=opts$message, complete=ESTIMATE$hess)
}
} else if(opts$SE.type == 'Fisher' && !(opts$method %in% c('MHRM', 'SEM', 'MIXED'))){
if(logPrior != 0 && opts$warn)
warning('Information matrix with the Fisher method does not
account for prior parameter distribution information')
ESTIMATE <- SE.Fisher(PrepList=PrepList, ESTIMATE=ESTIMATE, Theta=Theta, Data=Data,
constrain=constrain, Ls=Ls, full=opts$full,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND, warn=opts$warn,
omp_threads=opts$omp_threads)
}
ESTIMATE$cycles <- tmp$cycles
ESTIMATE$Prior <- tmp$Prior
ESTIMATE$Etable <- tmp$Etable
rm(tmp)
}
opts$times$end.time.SE <- proc.time()[3L]
opts$times$start.time.post <- proc.time()[3L]
cmods <- vector('list', Data$ngroups)
if(is.null(opts$theta_lim))
opts$theta_lim <- numeric(1)
lrPars <- ESTIMATE$lrPars
class(lrPars) <- 'S4'
for(g in seq_len(Data$ngroups)){
if(opts$method == 'MIXED' || opts$dentype == "discrete"){
F <- matrix(NA)
h2 <- numeric(1)
} else {
F <- Lambdas(ESTIMATE$pars[[g]], Names=colnames(data))
colnames(F) <- PrepList[[1L]]$factorNames
h2 <- rowSums(F^2)
}
cmods[[g]] <- new('SingleGroupClass', ParObjects=list(pars=ESTIMATE$pars[[g]], lrPars=lrPars,
random=ESTIMATE$random,
lr.random=ESTIMATE$lr.random),
Data = list(K=Data$K, nitems=Data$nitems),
Model = list(itemloc=PrepList[[1L]]$itemloc, nfact=nfact, constrain=constrain,
factorNames=PrepList[[1L]]$factorNames,
itemtype=PrepList[[1L]]$itemtype,
prodlist=PrepList[[1L]]$prodlist),
Options = list(method = 'MHRM', exploratory=PrepList[[1L]]$exploratory,
theta_lim=opts$theta_lim, dentype=opts$dentype),
Fit = list(G2=G2group[g], F=F, h2=h2),
Internals = list(Pl = rlist[[g]]$expected, CUSTOM.IND=CUSTOM.IND,
SLOW.IND=SLOW.IND))
if(opts$dentype %in% c("discrete", 'EH', 'EHW', 'Davidian', 'custom')){
cmods[[g]]@Model$Theta <- Theta
cmods[[g]]@Internals$Prior <- list(ESTIMATE$Prior[[g]])
}
}
#missing stats for MHRM
if(opts$method %in% c('MHRM', 'MIXED', 'SEM') &&
(!opts$logLik_if_converged || !(!ESTIMATE$converge && opts$logLik_if_converged))){
logLik <- G2 <- SElogLik <- 0
if(opts$draws > 0L){
if(opts$verbose) cat("\nCalculating log-likelihood...\n")
flush.console()
if(!opts$technical$parallel){
ncores <- .mirtClusterEnv$ncores
.mirtClusterEnv$ncores <- 1L
}
for(g in seq_len(Data$ngroups)){
cmods[[g]]@Data <- list(data=Data$data[Data$group == Data$groupName[g], ],
fulldata=Data$fulldata[[g]], tabdata=Data$tabdata,
Freq=list(Data$Freq[[g]]), K=Data$K)
cmods[[g]] <- calcLogLik(cmods[[g]], opts$draws, G2 = 'return',
lrPars=ESTIMATE$lrPars)
cmods[[g]]@Data <- list(K=Data$K, nitems=Data$nitems)
logLik <- logLik + cmods[[g]]@Fit$logLik
logPrior <- logPrior + cmods[[g]]@Fit$logPrior
SElogLik <- SElogLik + cmods[[g]]@Fit$SElogLik
G2 <- G2 + cmods[[g]]@Fit$G2
Pl[[g]] <- cmods[[g]]@Internals$Pl
}
if(!opts$technical$parallel)
.mirtClusterEnv$ncores <- ncores
}
}
####post estimation stats
if(opts$Moptim %in% c('solnp', 'nloptr')){
if(!is.null(opts$solnp_args$eqfun))
df <- df + length(opts$solnp_args$eqfun(ESTIMATE$shortpars, list()))
if(!is.null(opts$solnp_args$eval_g_eq))
df <- df + length(opts$solnp_args$eval_g_eq(ESTIMATE$shortpars, list()))
}
r <- rr
N <- sum(r)
tmp <- dfsubtr
AIC <- (-2) * logLik + 2 * tmp
BIC <- (-2) * logLik + tmp*log(N)
SABIC <- (-2) * logLik + tmp*log((N+2)/24)
HQ <- (-2) * logLik + 2*tmp*log(log(N))
p.G2 <- 1 - pchisq(G2,df)
RMSEA.G2 <- rmsea(X2=G2, df=df, N=N)
null.mod <- unclass(new('SingleGroupClass'))
TLI.G2 <- CFI.G2 <- NaN
if(opts$calcNull && length(r) * 3L < prod(Data$K) && opts$warn)
warning(c('Full table of responses is very sparse. ',
'Goodness-of-fit statistics may be very inaccurate'), call.=FALSE)
if(!opts$NULL.MODEL && opts$method != 'MIXED' && opts$calcNull && nmissingtabdata == 0L){
null.mod <- try(unclass(computeNullModel(data=data, key=key,
group=if(length(pars) > 1L) group else NULL)))
if(is(null.mod, 'try-error')){
if(opts$warn)
warning('Null model calculation did not converge.')
null.mod <- unclass(new('SingleGroupClass'))
} else if(!is.nan(G2)) {
TLI.G2 <- tli(X2=G2, X2.null=null.mod@Fit$G2, df=df, df.null=null.mod@Fit$df)
CFI.G2 <- cfi(X2=G2, X2.null=null.mod@Fit$G2, df=df, df.null=null.mod@Fit$df)
}
}
if(nmissingtabdata > 0L)
p.G2 <- RMSEA.G2 <- G2 <- TLI.G2 <- CFI.G2 <- NaN
if(is.null(parprior)) parprior <- list()
if(is.null(opts$quadpts)) opts$quadpts <- NaN
opts$times$end.time.post <- proc.time()[3L]
time <- opts$times
opts$times <- NULL
# set macro objects
Options <- opts
Options$exploratory <- PrepList[[1L]]$exploratory
Fit <- list(G2=G2, p=p.G2, TLI=TLI.G2, CFI=CFI.G2, RMSEA=RMSEA.G2, df=df,
AIC=AIC, BIC=BIC, SABIC=SABIC, HQ=HQ, logLik=logLik,
logPrior=logPrior, SElogLik=SElogLik, F=F, h2=h2)
pis <- if(opts$dentype == 'mixture')
ExtractMixtures(lapply(cmods, function(x) x@ParObjects$pars)) else NULL
Model <- list(model=oldmodel, factorNames=PrepList[[1L]]$factorNames, itemtype=PrepList[[1L]]$itemtype,
itemloc=PrepList[[1L]]$itemloc, nfact=nfact, pis=pis,
Theta=Theta, constrain=constrain, parprior=parprior, nest=as.integer(dfsubtr),
invariance=invariance, lrPars=lrPars, formulas=attr(mixed.design, 'formula'),
prodlist=PrepList[[1L]]$prodlist, nestpars=nestpars)
if(!is.null(opts$technical$Etable)){
Model$Etable <- ESTIMATE$rlist
Model$Etable$Theta <- Theta
}
Data$covdata <- if(length(lrPars)) lrPars@df else attr(mixed.design, 'covdata')
Data$itemdesign <- attr(mixed.design, 'itemdesign')
ParObjects <- list(pars=cmods, lrPars=lrPars, random=ESTIMATE$random, lr.random=ESTIMATE$lr.random)
OptimInfo <- list(iter=ESTIMATE$cycles, converged=ESTIMATE$converge, cand.t.var=ESTIMATE$cand.t.var,
condnum=NA, secondordertest=NA, SEMconv=SEMconv, aveAR=ESTIMATE$aveAR)
vcov <- matrix(NA, 1, 1)
if(Options$SE){
information <- ESTIMATE$info
if(!is.null(opts$technical$infoAsVcov)){
vcov <- information
} else {
if(!ESTIMATE$fail_invert_info){
isna <- is.na(diag(information))
info <- information[!isna, !isna]
vcov <- matrix(NA, ncol(information), ncol(information))
rownames(vcov) <- colnames(vcov) <- colnames(information)
vcov2 <- try(solve(info), silent=TRUE)
vcov[!isna, !isna] <- vcov2
if(!is(vcov2, 'try-error')){
OptimInfo$condnum <- kappa(info, exact=TRUE)
OptimInfo$secondordertest <- secondOrderTest(info)
} else OptimInfo$secondordertest <- FALSE
} else {
OptimInfo$secondordertest <- FALSE
}
}
} else OptimInfo$secondordertest <- NA
Internals <- list(collectLL=ESTIMATE$collectLL, Prior=ESTIMATE$Prior, Pl=Pl,
shortpars=as.numeric(ESTIMATE$shortpars), key=key,
bfactor=list(), CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND,
survey.weights=survey.weights, theta_lim = opts$theta_lim,
customGroup=customGroup, customItems=customItems)
if(opts$method == 'EM'){
tmp <- lapply(ESTIMATE$Etable, function(tab)
data.frame(Theta, posterior=rowSums(tab$r1)))
if(length(tmp)){
names(tmp) <- Data$groupNames
Internals$thetaPosterior <- tmp
}
}
if(opts$storeEtable)
Internals$Etable <- ESTIMATE$Etable
if(opts$storeEMhistory)
Internals$EMhistory <- ESTIMATE$EMhistory
if(opts$method == 'SEM') Options$TOL <- NA
if(opts$odentype == "discrete"){
Fit$F <- Fit$h2 <- NULL
mod <- new('DiscreteClass',
Data=Data,
Options=Options,
Fit=Fit,
Model=Model,
ParObjects=ParObjects,
OptimInfo=OptimInfo,
Internals=Internals,
vcov=vcov)
} else {
if(Data$ngroups == 1L){
ParObjects$pars <- cmods[[1L]]@ParObjects$pars
Internals$Pl <- Internals$Pl[[1L]]
if(opts$method == 'MIXED'){
Fit$p <- NaN
mod <- new('MixedClass',
Data=Data,
Options=Options,
Fit=Fit,
Model=Model,
ParObjects=ParObjects,
OptimInfo=OptimInfo,
Internals=Internals,
vcov=vcov)
} else {
if(Options$exploratory){
FF <- F %*% t(F)
V <- eigen(FF)$vector[ ,1L:nfact]
L <- eigen(FF)$values[1L:nfact]
if (nfact == 1L) F <- as.matrix(V * sqrt(L))
else F <- V %*% sqrt(diag(L))
if (sum(F[ ,1L] < 0)) F <- (-1) * F
colnames(F) <- paste("F", 1L:ncol(F), sep="")
h2 <- rowSums(F^2)
} else {
if(opts$method == 'EM')
Internals$bfactor <- list(prior=ESTIMATE$prior,
Priorbetween=ESTIMATE$Priorbetween,
sitems=ESTIMATE$sitems, specific=specific)
}
mod <- new('SingleGroupClass',
Data=Data,
Options=Options,
Fit=Fit,
Model=Model,
ParObjects=ParObjects,
OptimInfo=OptimInfo,
Internals=Internals,
vcov=vcov)
}
} else {
if(opts$method == 'EM')
Internals$bfactor <- list(prior=ESTIMATE$prior,
Priorbetween=ESTIMATE$Priorbetween,
sitems=ESTIMATE$sitems, specific=specific)
cls <- ifelse(opts$odentype == 'mixture', 'MixtureClass', 'MultipleGroupClass')
mod <- new(cls,
Data=Data,
Options=Options,
Fit=Fit,
Model=Model,
ParObjects=ParObjects,
OptimInfo=OptimInfo,
Internals=Internals,
vcov=vcov)
}
}
mod@time <- c("TOTAL:" = as.numeric(proc.time()[3L] - time$start.time),
Data = as.numeric(time$end.time.Data - time$start.time.Data),
ESTIMATE$time,
SE = as.numeric(time$end.time.SE - time$start.time.SE),
Post = as.numeric(time$end.time.post - time$start.time.post))
return(mod)
}
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 ESTIMATION
ESTIMATION <- function(data, model, group, itemtype = NULL, guess = 0, upper = 1,
invariance = '', pars = NULL, constrain = NULL, key = NULL,
parprior = NULL, mixed.design = NULL, customItems = NULL,
customItemsData = NULL, customGroup = NULL,
GenRandomPars = FALSE, large = FALSE,
survey.weights = NULL, latent.regression = NULL,
gpcm_mats=list(), spline_args=list(), monopoly.k=1,
control = list(), ...)
{
start.time <- proc.time()[3L]
dots <- list(...)
if(!is.null(itemtype))
itemtype <- ifelse(itemtype == 'grsm', 'grsmIRT', itemtype)
if(missing(data)) missingMsg('data')
if(length(unique(colnames(data))) != ncol(data))
stop('items must have unique names in data input', call.=FALSE)
if(missing(model)) missingMsg('model')
if(missing(group)) missingMsg('group')
if(!(is.factor(group) || is.character(group)) || length(group) != nrow(data))
stop('group input provided is not valid', call.=FALSE)
if(is.character(large) && large == 'return'){
Data <- opts <- list()
opts$zeroExtreme <- FALSE
opts$dentype <- 'default'
if(!is.null(dots$technical$zeroExtreme)) opts$zeroExtreme <- dots$technical$zeroExtreme
data <- as.matrix(data)
if(is.numeric(data))
data <- matrix(as.integer(data), nrow(data), dimnames=list(rownames(data), colnames(data)))
rownames(data) <- 1L:nrow(data)
if(is.null(colnames(data)))
colnames(data) <- paste0('Item.', 1L:ncol(data))
Data$data <- data
Data$group <- factor(group)
Data$groupNames <- levels(Data$group)
Data$ngroups <- length(Data$groupNames)
Data$nitems <- ncol(data)
K <- apply(Data$data, 2L, function(x) length(unique(na.omit(x))))
tmp <- list(...)
if(!is.null(tmp$technical$customK)) K <- tmp$technical$customK
itemloc <- c(1L, cumsum(K) + 1L)
Names <- NULL
for(i in 1L:length(K))
Names <- c(Names, paste("Item.",i,"_",1L:K[i],sep=""))
PrepListFull <- list(K=K, itemloc=itemloc,
Names=Names, itemnames=colnames(Data$data))
} else {
if(missing(data) || is.null(nrow(data)))
stop('data argument is required', call.=FALSE)
if(missing(model) || !is(model, 'numeric') && !is(model, 'mirt.model') &&
!is(model, 'character'))
stop('model argument (numeric, character, or from mirt.model() function) is required',
call.=FALSE)
if(!(is.factor(group) || is.character(group)) || length(group) != nrow(data))
stop('group input provided is not valid', call.=FALSE)
if(is.matrix(itemtype)){
itemtypefull <- itemtype
itemtype <- itemtypefull[1L,]
} else itemtypefull <- NULL
if(!is.null(itemtype))
stopifnot(is(itemtype, 'character'))
if(!is.null(constrain))
stopifnot(is(constrain, 'list'))
if(!is.null(parprior))
stopifnot(is(parprior, 'list'))
if(!is.null(customItems))
stopifnot(is(customItems, 'list'))
if(!is.null(customItemsData))
stopifnot(is(customItemsData, 'list') || length(customItemsData) == ncol(data))
if(!is.null(customGroup))
stopifnot(is(customGroup, 'GroupPars') || is.list(customGroup))
stopifnot(is(invariance, 'character'))
stopifnot(is(GenRandomPars, 'logical'))
stopifnot(is(large, 'logical') || is(large, 'list') || is(large, 'character'))
opts <- makeopts(GenRandomPars=GenRandomPars,
hasCustomGroup=!is.null(customGroup), ...)
if(opts$Moptim == 'NR'){
if(is.null(control$tol)) control$tol <- opts$TOL/1000
if(is.null(control$maxit)) control$maxit <- 50L
}
if(opts$Moptim == 'nlminb')
if(is.null(control$rel.tol)) control$rel.tol <- opts$TOL/100
if(opts$dentype == "discrete" && is.null(customGroup)){
den <- Theta_discrete_den
par <- if(!is.null(dots$structure)){
if(is.null(opts$technical$customTheta))
stop('customTheta must be defined when using structure input', call.=FALSE)
opts$structure <- model.matrix(dots$structure, as.data.frame(opts$technical$customTheta))
tmpnfact <- ncol(opts$technical$customTheta)
numeric(ncol(opts$structure))
} else if(is.null(opts$technical$customTheta)){
tmpnfact <- model
numeric(model-1L)
} else {
tmpnfact <- ncol(opts$technical$customTheta)
numeric(nrow(opts$technical$customTheta) - 1L)
}
if(length(par)){
names(par) <- paste0('c', 1:length(par))
est <- rep(TRUE, length(par))
} else {
par <- c(c = 0)
est <- FALSE
}
customGroup <- createGroup(par=par, est=est, den=den, nfact=tmpnfact,
gen=function(object) rnorm(length(object@par), 0, 1/2))
customGroup@itemclass <- -1L
rm(tmpnfact)
}
if(any(itemtype == 'ULL')){
opts$dentype <- 'custom'
den <- function(obj, Theta){
par <- obj@par
dlnorm(Theta, meanlog = par[1L], sdlog = par[2L])
}
opts$theta_lim <- c(.01, opts$theta_lim[2L]^2)
par <- c(meanlog=0, sdlog=1)
est <- c(FALSE, FALSE)
customGroup <- createGroup(par=par, est=est, den=den, nfact=1L,
gen=function(object) rnorm(length(object@par), 0, 1/2))
}
if(!is.null(survey.weights)){
stopifnot(opts$method %in% c('EM', 'QMCEM', 'MCEM'))
stopifnot(length(survey.weights) == nrow(data))
}
if(any(is.na(group))){
if(opts$message)
message('NA values in group removed, along with associated rows in data')
data <- data[!is.na(group), , drop=FALSE]
group <- group[!is.na(group)]
}
if(!is.null(customItems) || any(itemtype %in% Experimental_itemtypes()))
opts$calcNull <- FALSE
opts$times <- list(start.time=start.time)
# on exit, reset the seed to override internal
if(opts$method == 'MHRM' || opts$method == 'MIXED' || opts$method == 'SEM' &&
opts$plausible.draws == 0L)
on.exit(set.seed((as.numeric(Sys.time()) - floor(as.numeric(Sys.time()))) * 1e8))
#change itemtypes if NULL.MODEL
if(opts$NULL.MODEL){
constrain <- NULL
if(!is.null(itemtype)){
itemtype[itemtype == 'grsm'] <- 'graded'
itemtype[itemtype == 'rsm'] <- 'gpcm'
itemtype[itemtype == '3PL' | itemtype == '3PLu' | itemtype == '4PL'] <- '2PL'
itemtype[itemtype == '3PLNRM' | itemtype == '3PLuNRM' | itemtype == '4PLNRM'] <- '2PLNRM'
itemtype[itemtype == 'spline'] <- '2PL'
}
}
if(!is.null(itemtype)){
if(any(itemtype == 'spline') && !(opts$method %in% c('EM', 'QMCEM', 'MCEM')))
stop('spline itemtype only supported for EM algorithm', call.=FALSE)
}
if(length(group) != nrow(data))
stop('length of group not equal to number of rows in data.', call.=FALSE)
if(any(is.na(group)))
stop('Unknown group memberships cannot be estimated. Please remove the NA values in group
and the associated rows in the data input.', call.=FALSE)
#####
### uncomment when building html examples
# opts$verbose <- FALSE
####
opts$times$start.time.Data <- proc.time()[3L]
Data <- list()
data <- as.matrix(data)
if(!all(apply(data, 2L, class) %in% c('integer', 'numeric')))
stop('Input data must be integer or numeric values only', call.=FALSE)
if(is.numeric(data))
data <- matrix(as.integer(data), nrow(data),
dimnames=list(rownames(data), colnames(data)))
rownames(data) <- 1L:nrow(data)
if(is.null(colnames(data)))
colnames(data) <- paste0('Item.', 1L:ncol(data))
if(nrow(data) > 1L && is.null(opts$technical$customK)){
if(!is.null(key)){
key <- sapply(1L:length(key), function(ind, data, key){
if(is.na(key[ind])) return(NA)
s <- sort(unique(data[,ind]))
se <- min(s, na.rm = TRUE):(length(s) + min(s) - 1L)
ret <- se[s == key[ind]]
if(!length(ret)) ret <- NA
ret
}, data=data, key=key)
}
data <- remap.distance(data, message = opts$message)
}
Data$rowID <- 1L:nrow(data)
Data$completely_missing <- which(rowSums(is.na(data)) == ncol(data))
if(length(Data$completely_missing)){
if(opts$warn)
warning('data contains response patterns with only NAs', call.=FALSE)
pick <- rowSums(is.na(data)) != ncol(data)
Data$rowID <- unname(which(pick))
data <- subset(data, pick)
group <- subset(group, pick)
if(!is.null(latent.regression) || !is.null(mixed.design))
covdata <- covdata[-Data$completely_missing, , drop=FALSE]
}
Data$data <- data
if(is.null(opts$grsm.block)) Data$grsm.block <- rep(1L, ncol(data))
else Data$grsm.block <- opts$grsm.block
if(is.null(opts$rsm.block)) Data$rsm.block <- rep(1L, ncol(data))
else Data$rsm.block <- opts$rsm.block
Data$group <- factor(group)
Data$groupNames <- levels(Data$group)
if(any(grepl('-', Data$groupNames)))
stop('Group names cannot contain a dash (-) character', call.=FALSE)
Data$ngroups <- if(!is.null(opts$ngroups)) opts$ngroups else length(Data$groupNames)
if(!is.null(itemtypefull))
if(Data$ngroups != nrow(itemtypefull))
stop("number of groups does not match number of rows in itemtype", call.=FALSE)
Data$nitems <- ncol(Data$data)
Data$N <- nrow(Data$data)
Data$mins <- suppressWarnings(apply(data, 2L, min, na.rm=TRUE))
Data$mins[!is.finite(Data$mins)] <- 0L
if(!is.null(opts$technical$mins)) Data$mins <- opts$technical$mins
if(is.character(model)){
tmp <- any(sapply(colnames(data), grepl, x=model))
model <- mirt.model(model, itemnames = if(tmp) colnames(data) else NULL)
}
oldmodel <- model
PrepList <- vector('list', Data$ngroups)
tmp <- 1L:Data$ngroups
if(opts$dentype == "mixture")
Data$groupNames <- paste0("MIXTURE_", seq_len(Data$ngroups))
names(PrepList) <- Data$groupNames
model <- buildModelSyntax(model, J=Data$nitems, groupNames=Data$groupNames,
itemtype=itemtype)
Data$model <- model
if(!is.null(customGroup)){
if(Data$ngroups == 1L) customGroup <- list(customGroup)
else {
if(Data$ngroups != length(customGroup) && length(customGroup) == 1L){
tmplst <- vector('list', Data$ngroups)
for(g in seq_len(Data$ngroups))
tmplst[[g]] <- customGroup
customGroup <- tmplst
rm(tmplst)
}
}
}
if(!is.null(dots$PrepList)) {
PrepListFull <- PrepList[[1L]] <- dots$PrepList
} else {
PrepListFull <- PrepList[[1L]] <-
PrepData(data=Data$data, model=Data$model, itemtype=itemtype, guess=guess,
upper=upper, parprior=parprior, verbose=opts$verbose,
technical=opts$technical, parnumber=1L, BFACTOR=opts$dentype == 'bfactor',
grsm.block=Data$grsm.block, rsm.block=Data$rsm.block,
mixed.design=mixed.design, customItems=customItems,
customItemsData=customItemsData,
customGroup=customGroup[[1L]], spline_args=spline_args, monopoly.k=monopoly.k,
fulldata=opts$PrepList[[1L]]$fulldata, key=key, opts=opts,
gpcm_mats=gpcm_mats, internal_constraints=opts$internal_constraints,
dcIRT_nphi=opts$dcIRT_nphi, dentype=opts$dentype, item.Q=opts$item.Q)
if(!is.null(dots$Return_PrepList)) return(PrepListFull)
if(!is.null(itemtypefull)){
for(g in 2L:nrow(itemtypefull)){
PrepList[[g]] <-
PrepData(data=Data$data, model=Data$model, itemtype=itemtypefull[g,], guess=guess,
upper=upper, parprior=parprior, verbose=opts$verbose,
technical=opts$technical, parnumber=1L, BFACTOR=opts$dentype == 'bfactor',
grsm.block=Data$grsm.block, rsm.block=Data$rsm.block,
mixed.design=mixed.design, customItems=customItems,
customItemsData=customItemsData,
customGroup=customGroup[[1L]], spline_args=spline_args, monopoly.k=monopoly.k,
fulldata=opts$PrepList[[1L]]$fulldata, key=key, opts=opts,
gpcm_mats=gpcm_mats, internal_constraints=opts$internal_constraints,
dcIRT_nphi=opts$dcIRT_nphi, dentype=opts$dentype, item.Q=opts$item.Q)
}
}
}
if(!is.null(opts$structure)){
PrepList[[1L]]$pars[[Data$nitems+1L]]@structure <- opts$structure
PrepList[[1L]]$pars[[Data$nitems+1L]]@parnames <-
names(PrepList[[1L]]$pars[[Data$nitems+1L]]@est) <- colnames(opts$structure)
PrepList[[1L]]$pars[[Data$nitems+1L]]@itemclass <- -999L
}
parnumber <- max(PrepList[[1L]]$pars[[Data$nitems+1L]]@parnum) + 1L
attr(PrepListFull$pars, 'nclasspars') <- attr(PrepList[[1L]]$pars, 'nclasspars') <-
matrix(sapply(PrepListFull$pars, function(y) length(y@parnum)), nrow=1L)
for(g in seq_len(Data$ngroups)){
if(g != 1L){
if(is.null(itemtypefull))
PrepList[[g]] <- list(pars=PrepList[[1L]]$pars)
else attr(PrepList[[g]]$pars, 'nclasspars') <-
sapply(PrepList[[g]]$pars, function(y) length(y@parnum))
for(i in 1L:length(PrepList[[g]]$pars)){
PrepList[[g]]$pars[[i]]@parnum <- parnumber:(parnumber +
length(PrepList[[g]]$pars[[i]]@parnum) - 1L)
parnumber <- max(PrepList[[g]]$pars[[i]]@parnum) + 1L
}
if(!is.null(customGroup)){
customGroup[[g]]@parnum <- PrepList[[g]]$pars[[Data$nitems + 1L]]@parnum
PrepList[[g]]$pars[[Data$nitems + 1L]] <- customGroup[[g]]
}
}
}
if(length(mixed.design$random) > 0L){
for(i in seq_len(length(mixed.design$random))){
mixed.design$random[[i]]@parnum <- parnumber:(parnumber - 1L +
length(mixed.design$random[[i]]@par))
parnumber <- max(mixed.design$random[[i]]@parnum) + 1L
}
}
if(!is.null(latent.regression)){
if(length(PrepListFull$prodlist))
stop('Polynomial combinations currently not supported when latent regression effects are used', call.=FALSE)
lrPars <- make.lrdesign(df=latent.regression$df, formula=latent.regression$formula,
factorNames=PrepListFull$factorNames, EM=latent.regression$EM,
TOL=opts$TOL)
lrPars@parnum <- parnumber:(parnumber - 1L + length(lrPars@par))
parnumber <- max(lrPars@parnum) + 1L
if(opts$dentype == 'discrete'){
tmp <- matrix(1L:length(lrPars@beta), nrow(lrPars@beta), ncol(lrPars@beta))
tmp2 <- tmp[1, ]
lrPars@est[tmp2[-length(tmp2)]] <- TRUE
lrPars@est[tmp[,ncol(tmp)]] <- FALSE
PrepList$all$pars[[ncol(data) + 1L]]@est <- FALSE
}
} else lrPars <- list()
if(length(latent.regression$lr.random) > 0L){
for(i in seq_len(length(latent.regression$lr.random))){
latent.regression$lr.random[[i]]@parnum <- parnumber:(parnumber - 1L +
length(latent.regression$lr.random[[i]]@par))
parnumber <- max(latent.regression$lr.random[[i]]@parnum) + 1L
}
}
}
if(!is.null(opts$PrepList)){
PrepList <- opts$PrepList
} else {
if(!is.list(large)){
tmpdata <- Data$data
if(opts$zeroExtreme){
n_is_na <- rowSums(is.na(tmpdata))
sums <- colSums(t(tmpdata) - Data$mins, na.rm=TRUE)
maxs <- apply(tmpdata, 2L, function(x) length(unique(na.omit(x)))) - 1L
if(is.null(survey.weights)) survey.weights <- rep(1, nrow(tmpdata))
survey.weights[sums == 0L | (sums + n_is_na) == sum(maxs)] <- 0
}
tmptabdata <- if(is.logical(large) && large){
maketabDataLarge(tmpdata=tmpdata, group=Data$group,
groupNames=if(opts$dentype != 'mixture') Data$groupNames else 'full',
nitem=Data$nitems, K=PrepListFull$K, itemloc=PrepListFull$itemloc,
Names=PrepListFull$Names, itemnames=PrepListFull$itemnames,
survey.weights=survey.weights)
} else {
maketabData(tmpdata=tmpdata, group=Data$group,
groupNames=if(opts$dentype != 'mixture') Data$groupNames else 'full',
nitem=Data$nitems, K=PrepListFull$K, itemloc=PrepListFull$itemloc,
Names=PrepListFull$Names, itemnames=PrepListFull$itemnames,
survey.weights=survey.weights)
}
if(is.character(large) && large == 'return'){
return(tmptabdata)
} else large <- tmptabdata
rm(tmpdata, tmptabdata)
}
Data$tabdatalong <- large$tabdata
Data$tabdata <- large$tabdata2
for(g in seq_len(Data$ngroups)){
if(opts$dentype == 'mixture'){
Data$fulldata[[g]] <- PrepListFull$fulldata
Data$Freq[[g]] <- large$Freq[[1L]]
} else {
Data$fulldata[[g]] <- PrepListFull$fulldata[Data$group == Data$groupNames[g],
, drop=FALSE]
Data$Freq[[g]] <- large$Freq[[g]]
}
}
}
if(opts$returnPrepList) return(PrepList)
if(opts$dentype == 'bfactor'){
#better start values
if((PrepList[[1L]]$nfact - attr(model, 'nspec')) == 1L){
nfact <- PrepListFull$nfact
for(g in seq_len(Data$ngroups)){
for(i in seq_len(Data$nitems)){
tmp <- PrepList[[g]]$pars[[i]]@par[1L:nfact]
tmp2 <- tmp[1L]
tmp[PrepList[[g]]$pars[[i]]@est[1L:nfact]] <-
tmp[PrepList[[g]]$pars[[i]]@est[1L:nfact]]/2
tmp[1L] <- tmp2
PrepList[[g]]$pars[[i]]@par[1L:nfact] <- tmp
}
}
}
}
PrepList <- UpdateParameters(PrepList, model, groupNames=Data$groupNames)
if(GenRandomPars){
for(g in seq_len(Data$ngroups))
for(i in seq_len(length(PrepList[[g]]$pars)))
PrepList[[g]]$pars[[i]] <- GenRandomPars(PrepList[[g]]$pars[[i]])
}
if(opts$dentype == "discrete"){
PrepList[[1L]]$exploratory <- FALSE
if(is.null(opts$technical$customTheta))
opts$technical$customTheta <- diag(PrepList[[1L]]$nfact)
}
RETURNVALUES <- FALSE
SUPPLIED_STARTS <- FALSE
if(!is.null(pars)){
if(is(pars, 'data.frame')){
SUPPLIED_STARTS <- TRUE
PrepList <- UpdatePrepList(PrepList, pars, random=mixed.design$random,
lrPars=lrPars, lr.random=latent.regression$lr.random,
MG = TRUE)
mixed.design$random <- attr(PrepList, 'random')
latent.regression$lr.random <- attr(PrepList, 'lr.random')
if(any(pars$class == 'lrPars')) lrPars <- update.lrPars(pars, lrPars)
attr(PrepList, 'random') <- NULL
attr(PrepList, 'lr.random') <- NULL
}
if(!is.null(attr(pars, 'values')) || (is.character(pars)))
RETURNVALUES <- TRUE
}
pars <- vector('list', Data$ngroups)
for(g in seq_len(Data$ngroups))
pars[[g]] <- PrepList[[g]]$pars
nitems <- Data$nitems
Data$K <- PrepList[[1L]]$K
nfact <- PrepList[[1L]]$pars[[nitems+1L]]@nfact
if(nfact != 1L && any(c('Rasch') %in% itemtype ) && PrepList[[1L]]$exploratory)
stop('Rasch itemtype is for confirmatory models only.', call.=FALSE)
nLambdas <- PrepList[[1L]]$pars[[1L]]@nfact
if(is.null(constrain)) constrain <- list()
nspec <- ifelse(!is.null(attr(model, 'nspec')), attr(model, 'nspec'), 1L)
#default MG uses configural model (independent groups but each identified)
if('free_means' %in% invariance ){ #Free factor means (means 0 for ref)
if(opts$dentype == 'bfactor'){
for(g in 2L:Data$ngroups)
pars[[g]][[nitems + 1L]]@est[1L:(nfact-nspec)] <- TRUE
} else {
for(g in 2L:Data$ngroups)
pars[[g]][[nitems + 1L]]@est[1L:nfact] <- TRUE
}
}
dummymat <- matrix(FALSE, pars[[1L]][[nitems + 1L]]@nfact, pars[[1L]][[nitems + 1L]]@nfact)
if(any('free_means' %in% invariance)){ #Free means
if(all(sapply(PrepList[[1]]$pars, function(x) class(x)) %in%
c(ordinal_itemtypes(), 'GroupPars'))){
TS <- rowMeans(Data$data, na.rm=TRUE)
gmus <- tapply(TS, Data$group, mean, na.rm=TRUE)
gsd <- sd(TS[Data$group == Data$groupNames[1L]], na.rm=TRUE)
gmuscaled <- (gmus - gmus[Data$groupNames[1L]]) / gsd
if(all(is.finite(gmuscaled))){
for(i in 1L:Data$ngroups){
tmp <- pars[[i]][[Data$nitems+1L]]
tmp@par[tmp@est & grepl('MEAN_', names(tmp@est))] <- gmuscaled[i]
pars[[i]][[Data$nitems+1L]] <- tmp
}
}
}
}
if(any(c('free_var', 'free_vars') %in% invariance)){ #Free factor vars (vars 1 for ref)
if(opts$dentype == 'bfactor'){
tmp <- dummymat[1L:(nfact-nspec),1L:(nfact-nspec), drop=FALSE]
diag(tmp) <- TRUE
dummymat[1L:(nfact-nspec),1L:(nfact-nspec)] <- tmp
} else {
diag(dummymat) <- TRUE
}
tmp <- dummymat[lower.tri(dummymat, TRUE)]
if(opts$dentype %in% c('Davidian', 'EHW')){
for(g in 2L:Data$ngroups)
pars[[g]][[nitems + 1L]]@est[2L] <- TRUE
} else {
for(g in 2L:Data$ngroups){
pars[[g]][[nitems + 1L]]@est <- pars[[g]][[nitems + 1L]]@est |
c(pars[[g]][[nitems + 1L]]@est[1L:pars[[g]][[nitems + 1L]]@nfact], tmp)
names(pars[[g]][[nitems + 1L]]@est) <- names(pars[[g]][[nitems + 1L]]@par)
pars[[g]][[nitems + 1L]]@parnames <- names(pars[[g]][[nitems + 1L]]@est)
}
}
}
if(opts$dentype == 'mixture' && !SUPPLIED_STARTS){
tmp <- length(pars[[1L]][[nitems + 1L]]@par)
pars[[1L]][[nitems + 1L]]@est[tmp] <- FALSE
for(g in 1L:Data$ngroups)
pars[[g]][[nitems + 1L]]@par[tmp] <- g - 1
names(PrepList) <- Data$groupNames
}
constrain <- UpdateConstrain(pars=pars, constrain=constrain, invariance=invariance, nfact=Data$nfact,
nLambdas=nLambdas, J=nitems, ngroups=Data$ngroups, PrepList=PrepList,
method=opts$method, itemnames=PrepList[[1L]]$itemnames, model=model,
groupNames=Data$groupNames)
pars <- resetPriorConstrain(pars=pars, constrain=constrain,
nconstrain=opts$technical$nconstrain)
if(RETURNVALUES){
for(g in seq_len(Data$ngroups))
PrepList[[g]]$pars <- pars[[g]]
return(ReturnPars(PrepList, PrepList[[1L]]$itemnames, lr.random=latent.regression$lr.random,
random=mixed.design$random, lrPars=lrPars, MG = TRUE))
}
startlongpars <- c()
if(opts$NULL.MODEL){
constrain <- list()
for(g in seq_len(length(pars))){
for(i in seq_len(nitems)){
pars[[g]][[i]] <- set_null_model(pars[[g]][[i]])
}
}
}
rr <- 0L
for(g in seq_len(Data$ngroups)){
if(g > 1L && opts$dentype == 'mixture') break
r <- Data$Freq[[g]]
rr <- rr + r
}
df <- if(opts$dentype == 'mixture') prod(PrepList[[1L]]$K) - 1
else Data$ngroups * (prod(PrepList[[1L]]$K) - 1)
if(df > 1e10) df <- 1e10
nestpars <- nconstr <- 0L
for(g in seq_len(Data$ngroups))
for(i in seq_len(nitems+1L))
nestpars <- nestpars + sum(pars[[g]][[i]]@est)
if(!is.null(mixed.design$random)){
for(i in seq_len(length(mixed.design$random)))
nestpars <- nestpars + sum(mixed.design$random[[i]]@est)
}
if(length(lrPars))
nestpars <- nestpars + sum(lrPars@est)
if(!is.null(dots$structure)) nestpars <- nestpars - 1L
if(!is.null(opts$technical$nconstrain))
constrain <- c(constrain, opts$technical$nconstrain)
if(length(constrain) > 0L)
for(i in seq_len(length(constrain)))
nconstr <- nconstr + length(constrain[[i]]) - 1L
if(Data$ngroups > 1L && !length(constrain)){
if(opts$warn && any(invariance %in% c('free_means', 'free_var')))
warning('Multiple-group model may not be identified without providing anchor items',
call.=FALSE)
for(j in seq_len(Data$ngroups)){
if(opts$dentype != 'mixture')
tmp <- apply(subset(Data$data, Data$group == Data$groupNames[j]), 2L,
function(x) length(unique(na.omit(x)))) == Data$K
for(i in which(!tmp)){
if(any(PrepList[[j]]$pars[[i]]@est))
stop(paste0('Multiple Group model will not be identified without ',
'proper constraints (groups contain missing data patterns ',
'where item responses have been completely omitted or, alternatively, ',
'the number of categories within each group is not equal to the ',
'total number of categories)'),
call. = FALSE)
}
}
}
nmissingtabdata <- sum(is.na(rowSums(Data$tabdata)))
dfsubtr <- nestpars - nconstr
if(opts$dentype == 'EH') dfsubtr <- dfsubtr + (opts$quadpts - 1L) * Data$ngroups
else if(opts$dentype == 'EHW') dfsubtr <- dfsubtr + (opts$quadpts - 3L) * Data$ngroups
if(df <= dfsubtr)
stop('Too few degrees of freedom. There are only ', df, ' degrees of freedom but ',
dfsubtr, ' parameters were freely estimated.', call.=FALSE)
df <- df - dfsubtr
if(!is.null(customItems)){
for(g in seq_len(Data$ngroups))
PrepList[[g]]$exploratory <- FALSE
}
G2group <- numeric(Data$ngroups)
DERIV <- vector('list', Data$ngroups)
for(g in seq_len(Data$ngroups)){
DERIV[[g]] <- vector('list', Data$nitems)
for(i in seq_len(Data$nitems))
DERIV[[g]][[i]] <- selectMethod(Deriv, c(class(pars[[g]][[i]]), 'matrix'))
}
Ls <- makeLmats(pars, constrain, random = mixed.design$random,
lr.random=latent.regression$lr.random, lrPars=lrPars,
nconstrain=opts$technical$nconstrain)
CUSTOM.IND <- which(sapply(pars[[1L]], class) %in% Use_R_ProbTrace())
SLOW.IND <- which(sapply(pars[[1L]], class) %in% Use_R_Deriv())
if(pars[[1]][[length(pars[[1L]])]]@itemclass %in% c(-1L, -999L))
SLOW.IND <- c(SLOW.IND, length(pars[[1L]]))
if(opts$dentype != 'Gaussian' && opts$method %in% c('MHRM', 'MIXED', 'SEM'))
stop('Non-Gaussian densities not currently supported with MHRM algorithm')
#warnings
wmsg <- 'Lower and upper bound parameters (g and u) should use \'norm\' (i.e., logit) prior'
for(g in seq_len(length(pars))){
for(i in seq_len(length(pars[[1L]]))){
if(is(pars[[g]][[i]], 'dich')){
pt <- pars[[g]][[i]]@prior.type
if(!(pt[length(pt)-1L] %in% c(0L, 1L, 4L))) warning(wmsg, call.=FALSE)
if(!(pt[length(pt)] %in% c(0L, 1L, 4L))) warning(wmsg, call.=FALSE)
next
} else if(is(pars[[g]][[i]], 'partcomp')){
pt <- pars[[g]][[i]]@prior.type
if(!(pt[length(pt)-1L] %in% c(0L, 1L))) warning(wmsg, call.=FALSE)
if(!(pt[length(pt)] %in% c(0L, 1L))) warning(wmsg, call.=FALSE)
next
} else if(is(pars[[g]][[i]], 'nestlogit')){
pt <- pars[[g]][[i]]@prior.type
if(!(pt[nfact + 2L] %in% c(0L, 1L))) warning(wmsg, call.=FALSE)
if(!(pt[nfact + 3L] %in% c(0L, 1L))) warning(wmsg, call.=FALSE)
next
}
}
}
SEMconv <- NA
Data$wmiss <- if(length(lrPars)) with(Data, 1/rowMeans(!is.na(data)) / nitems)
else with(Data, 1/rowMeans(!is.na(tabdata)) / nitems)
opts$times$end.time.Data <- proc.time()[3L]
#EM estimation
opts$times$start.time.Estimate <- proc.time()[3L]
logLik <- G2 <- SElogLik <- NaN
logPrior <- 0
Pl <- vector('list', Data$ngroups)
if(opts$method %in% c('EM', 'BL', 'QMCEM', 'MCEM')){
logLik <- G2 <- SElogLik <- 0
if(length(lrPars)){
if(opts$SE && !(opts$SE.type %in% c('complete', 'forward', 'central', 'Richardson')))
stop('Information matrix method for latent regression estimates not supported',
call.=FALSE)
opts$full <- TRUE
} else opts$full <- FALSE
temp <- matrix(0L,nrow=nitems,ncol=nspec)
sitems <- matrix(0L, nrow=sum(PrepList[[1L]]$K), ncol=nspec)
specific <- NULL
if(opts$dentype == "discrete"){
theta <- 0
Theta <- opts$technical$customTheta
opts$quadpts <- nrow(Theta)
if(pars[[1L]][[1L]]@nfact != ncol(Theta))
stop("mirt.model definition does not have same number of traits/attributes as customTheta input", call.=FALSE)
} else {
if(is.null(opts$quadpts)){
tmp <- if(opts$dentype == 'bfactor') PrepList[[1L]]$nfact - attr(model, 'nspec') + 1L
else nfact
opts$quadpts <- select_quadpts(tmp)
}
if(opts$quadpts < 3 && opts$warn) warning('Should use more than 2 quadpts', call.=FALSE)
theta <- 1
if(!(opts$method %in% c('QMCEM', 'MCEM')))
theta <- as.matrix(seq(opts$theta_lim[1L], opts$theta_lim[2L],
length.out = opts$quadpts))
if(opts$dentype == 'bfactor'){
specific <- attr(oldmodel, 'specific')
specific[is.na(specific)] <- 1L
for(i in seq_len(nitems)) temp[i, specific[i]] <- 1L
ind <- 1L
for(i in seq_len(nitems)){
for(j in seq_len(PrepList[[1L]]$K[i])){
sitems[ind, ] <- temp[i, ]
ind <- ind + 1L
}
}
nfact2 <- PrepList[[1L]]$nfact - attr(model, 'nspec') + 1L
Theta <- thetaComb(theta, nfact2)
Theta <- cbind(Theta[,1L:(nfact2-1L),drop=FALSE],
matrix(Theta[,nfact2], nrow=nrow(Theta), ncol=ncol(sitems)))
} else {
if(opts$method %in% c('QMCEM', 'MCEM')){
Theta <- NULL
} else {
if(opts$quadpts^nfact <= opts$MAXQUAD){
if(is.null(opts$technical$customTheta))
Theta <- thetaComb(theta, nfact)
} else stop('Greater than ', opts$MAXQUAD, ' quadrature points.', call.=FALSE)
if(opts$message && nfact > 3L && !(opts$odentype %in% c('custom', 'discrete')))
message('EM quadrature for high dimensional models are better handled
\twith the \"QMCEM\" or \"MCEM\" method')
}
}
if(!is.null(opts$technical$customTheta)){
Theta <- opts$technical$customTheta
if(!is.matrix(Theta)) stop('customTheta input must be a matrix', call.=FALSE)
opts$quadpts <- nrow(Theta)
}
pars <- loadSplinePars(pars, Theta)
} #end Theta def
ESTIMATE <- EM.group(pars=pars, constrain=constrain, Ls=Ls, PrepList=PrepList, Data=Data,
list = list(NCYCLES=opts$NCYCLES, TOL=opts$TOL, MSTEPTOL=opts$MSTEPTOL,
nfactNames=PrepList[[1L]]$nfactNames, theta=theta,
itemloc=PrepList[[1L]]$itemloc, dentype=opts$dentype,
sitems=sitems, specific=specific, NULL.MODEL=opts$NULL.MODEL,
nfact=nfact, constrain=constrain, verbose=opts$verbose,
SE = opts$SE, SE.type=opts$SE.type, delta=opts$delta, quadpts=opts$quadpts,
accelerate=opts$accelerate, CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND,
customPriorFun=opts$customPriorFun, Moptim=opts$Moptim, warn=opts$warn,
message=opts$message, method=opts$method, full=opts$full,
lrPars=lrPars, SE=opts$SE && opts$SE.type == 'numerical', Etable=opts$Etable,
NULL.MODEL=opts$NULL.MODEL, PLCI=opts$PLCI, Norder=opts$Norder,
keep_vcov_PD=opts$keep_vcov_PD, symmetric=opts$technical$symmetric,
MCEM_draws=opts$MCEM_draws, omp_threads=opts$omp_threads),
Theta=Theta, DERIV=DERIV, solnp_args=opts$solnp_args, control=control,
nconstrain=opts$technical$nconstrain)
if(opts$method == 'MCEM')
opts$quadpts <- opts$MCEM_draws(ESTIMATE$cycles)
opts$Moptim <- ESTIMATE$Moptim
lrPars <- ESTIMATE$lrPars
startlongpars <- ESTIMATE$longpars
rlist <- ESTIMATE$rlist
logPrior <- ESTIMATE$logPrior
for(g in seq_len(Data$ngroups)){
if(g > 1L && opts$dentype == 'mixture') break
Pl[[g]] <- rlist[[g]]$expected
Pltmp <- Pl[[g]]
if(opts$full){
rg <- 1
G2group[g] <- NaN
} else {
rg <- Data$Freq[[g]]
Pltmp <- Pltmp[rg != 0]
rg <- rg[rg != 0]
Ng <- sum(rg)
G2group[g] <- 2 * sum(rg * log(rg/(Ng*Pltmp)))
}
G2 <- G2 + G2group[g]
logLik <- logLik + sum(rg*log(Pltmp))
}
} else if(opts$method %in% c('MHRM', 'SEM')){ #MHRM estimation
Theta <- matrix(0, Data$N, nitems)
if(opts$method == 'SEM') opts$NCYCLES <- NA
ESTIMATE <- MHRM.group(pars=pars, constrain=constrain, Ls=Ls, PrepList=PrepList, Data=Data,
list = list(NCYCLES=opts$NCYCLES, BURNIN=opts$BURNIN,
SEMCYCLES=opts$SEMCYCLES, gain=opts$gain,
KDRAWS=opts$KDRAWS, MHDRAWS=opts$MHDRAWS,
TOL=opts$TOL, SE=FALSE, SE.type = 'none',
nfactNames=PrepList[[1L]]$nfactNames,
itemloc=PrepList[[1L]]$itemloc,
nfact=nfact, constrain=constrain, verbose=opts$verbose,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND,
startlongpars=startlongpars,
cand.t.var=opts$technical$MHcand, warn=opts$warn,
message=opts$message, expl=PrepList[[1L]]$exploratory,
plausible.draws=opts$plausible.draws,
MSTEPTOL=opts$MSTEPTOL, Moptim=opts$Moptim,
keep_vcov_PD=opts$keep_vcov_PD),
DERIV=DERIV, solnp_args=opts$solnp_args, control=control)
if(opts$plausible.draws != 0) return(ESTIMATE)
if(opts$SE && (ESTIMATE$converge || !opts$info_if_converged)){
if(opts$verbose)
cat('\nCalculating information matrix...\n')
tmp <- MHRM.group(pars=ESTIMATE$pars, constrain=constrain, Ls=Ls, PrepList=PrepList, Data=Data,
list = list(NCYCLES=opts$MHRM_SE_draws, BURNIN=1L,
SEMCYCLES=opts$SEMCYCLES, gain=opts$gain,
KDRAWS=opts$KDRAWS, MHDRAWS=opts$MHDRAWS,
TOL=opts$SEtol, SE=TRUE, SE.type=opts$SE.type,
nfactNames=PrepList[[1L]]$nfactNames,
itemloc=PrepList[[1L]]$itemloc,
nfact=nfact, constrain=constrain, verbose=FALSE,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND,
startlongpars=ESTIMATE$longpars, plausible.draws=0L,
cand.t.var=opts$technical$MHcand, warn=opts$warn,
message=opts$message, expl=PrepList[[1L]]$exploratory,
MSTEPTOL=opts$MSTEPTOL, Moptim='NR1',
keep_vcov_PD=opts$keep_vcov_PD),
DERIV=DERIV, solnp_args=opts$solnp_args, control=control)
ESTIMATE$pars <- tmp$pars
ESTIMATE$info <- tmp$info
ESTIMATE$fail_invert_info <- tmp$fail_invert_info
ESTIMATE$time <- c(ESTIMATE$time, SE=sum(tmp$time))
rm(tmp)
}
rlist <- vector('list', Data$ngroups)
for(g in seq_len(Data$ngroups))
rlist[[g]]$expected = numeric(1L)
} else if(opts$method == 'MIXED'){
if(is.null(opts$technical$RANDSTART)) opts$technical$RANDSTART <- 100L
if(is.null(opts$technical$BURNIN) && length(mixed.design$random)) opts$BURNIN <- 200L
Theta <- matrix(0, Data$N, nitems)
ESTIMATE <- MHRM.group(pars=pars, constrain=constrain, Ls=Ls,
PrepList=PrepList, random=mixed.design$random, Data=Data,
lrPars=lrPars, lr.random=latent.regression$lr.random,
list = list(NCYCLES=opts$NCYCLES, BURNIN=opts$BURNIN,
SEMCYCLES=opts$SEMCYCLES, gain=opts$gain,
KDRAWS=opts$KDRAWS, MHDRAWS=opts$MHDRAWS,
TOL=opts$TOL, SE.type = 'none',
nfactNames=PrepList[[1L]]$nfactNames,
itemloc=PrepList[[1L]]$itemloc,
nfact=nfact, constrain=constrain, verbose=opts$verbose,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND,
startlongpars=startlongpars, SE=FALSE,
cand.t.var=opts$technical$MHcand, warn=opts$warn,
message=opts$message, expl=FALSE, plausible.draws=0L,
RANDSTART=opts$technical$RANDSTART,
MSTEPTOL=opts$MSTEPTOL, Moptim=opts$Moptim,
keep_vcov_PD=opts$keep_vcov_PD),
DERIV=DERIV, solnp_args=opts$solnp_args, control=control)
if(opts$SE && (ESTIMATE$converge || !opts$info_if_converged)){
if(opts$verbose)
cat('\nCalculating information matrix...\n')
tmp <- MHRM.group(pars=ESTIMATE$pars, constrain=constrain, Ls=Ls,
PrepList=PrepList, random=mixed.design$random, Data=Data,
lrPars=ESTIMATE$lrPars, lr.random=latent.regression$lr.random,
list = list(NCYCLES=opts$MHRM_SE_draws, BURNIN=1L,
SEMCYCLES=opts$SEMCYCLES, gain=opts$gain,
KDRAWS=opts$KDRAWS, MHDRAWS=opts$MHDRAWS,
TOL=opts$SEtol, SE=TRUE, SE.type=opts$SE.type,
nfactNames=PrepList[[1L]]$nfactNames,
itemloc=PrepList[[1L]]$itemloc,
nfact=nfact, constrain=constrain, verbose=FALSE,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND,
startlongpars=ESTIMATE$longpars, plausible.draws=0L,
cand.t.var=opts$technical$MHcand, warn=opts$warn,
message=opts$message, expl=FALSE,
RANDSTART=1L,
MSTEPTOL=opts$MSTEPTOL, Moptim='NR1',
keep_vcov_PD=opts$keep_vcov_PD),
DERIV=DERIV, solnp_args=opts$solnp_args, control=control)
ESTIMATE$pars <- tmp$pars
ESTIMATE$random <- tmp$random
ESTIMATE$lrPars <- tmp$lrPars
ESTIMATE$lr.random <- tmp$lr.random
ESTIMATE$info <- tmp$info
ESTIMATE$fail_invert_info <- tmp$fail_invert_info
ESTIMATE$time <- c(ESTIMATE$time, SE=sum(tmp$time))
rm(tmp)
}
rlist <- vector('list', Data$ngroups)
for(g in seq_len(Data$ngroups))
rlist[[g]]$expected = numeric(1L)
}
for(g in seq_len(length(pars))){
for(i in seq_len(length(pars[[1L]]))){
if(is(pars[[g]][[i]], 'dich')){
tmp <- ESTIMATE$pars[[g]][[i]]@par
nms <- ESTIMATE$pars[[g]][[i]]@parnames
if(tmp[nms == 'g'] > tmp[nms == 'u']){
if(opts$warn)
warning('g parameter greater than u detected. Model did not converge', call.=FALSE)
ESTIMATE$converge <- FALSE
}
}
}
}
opts$times$end.time.Estimate <- proc.time()[3L]
if(opts$logLik_if_converged && !ESTIMATE$converge) opts$draws <- 0
opts$times$start.time.SE <- proc.time()[3L]
if(!opts$NULL.MODEL && opts$SE){
tmp <- ESTIMATE
if(opts$verbose && !(opts$method %in% c('MHRM', 'MIXED', 'SEM')))
cat('\n\nCalculating information matrix...\n')
if(opts$SE.type %in% c('complete', 'Oakes') && opts$method %in% c('EM', 'QMCEM')){
opts$times$start.time.SE <- ESTIMATE$start.time.SE
ESTIMATE <- loadESTIMATEinfo(info=-ESTIMATE$hess, ESTIMATE=ESTIMATE, constrain=constrain,
warn=opts$warn)
} else if(opts$SE.type == 'SEM' && opts$method == 'EM'){
collectLL <- as.numeric(ESTIMATE$collectLL)
collectLL <- exp(c(NA, collectLL) - c(collectLL, NA))
from <- suppressWarnings(max(which(collectLL <= opts$SEM_from)))
if(from < 1L) from <- 1L
to <- min(which(collectLL >= opts$SEM_to))
dontrun <- FALSE
if(from == to){
if(opts$warn)
warning('SEM window is too small to compute information matrix.
Consider changing the starting values', call.=FALSE)
dontrun <- TRUE
}
lengthsplit <- do.call(c, lapply(strsplit(names(ESTIMATE$correct), 'COV_'), length))
lengthsplit <- lengthsplit + do.call(c, lapply(strsplit(names(ESTIMATE$correct), 'MEAN_'), length))
is.latent <- lengthsplit > 2L
if(!dontrun){
if(ESTIMATE$cycles <= 10L)
if(opts$warn)
warning('Very few EM cycles performed. Consider decreasing TOL further to
increase EM iteration count or starting farther away from ML estimates by
passing the \'GenRandomPars = TRUE\' argument')
estmat <- matrix(FALSE, length(ESTIMATE$correction), length(ESTIMATE$correction))
DM <- estmat + 0
diag(estmat) <- TRUE
if(!opts$technical$parallel){
ncores <- .mirtClusterEnv$ncores
.mirtClusterEnv$ncores <- 1L
}
DM <- myLapply(1L:ncol(estmat), FUN=SE.SEM, estmat=estmat, pars=ESTIMATE$pars, constrain=constrain, Data=Data,
list = list(NCYCLES=opts$NCYCLES, TOL=opts$SEtol, MSTEPTOL=opts$MSTEPTOL,
nfactNames=PrepList[[1L]]$nfactNames, theta=theta,
itemloc=PrepList[[1L]]$itemloc, keep_vcov_PD=opts$keep_vcov_PD,
sitems=sitems, specific=specific, NULL.MODEL=opts$NULL.MODEL,
nfact=nfact, constrain=constrain, verbose=opts$verbose,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND, Moptim=ESTIMATE$Moptim,
EHPrior=ESTIMATE$Prior, warn=opts$warn, dentype=opts$dentype,
message=opts$message, full=opts$full, lrPars=lrPars, method=opts$method),
Theta=Theta, theta=theta, ESTIMATE=ESTIMATE, from=from, to=to,
DERIV=DERIV, is.latent=is.latent, Ls=Ls, PrepList=PrepList,
solnp_args=opts$solnp_args, control=control)
SEMconv <- sapply(DM, function(x) all(attr(x, 'converged')))
if(!all(SEMconv)){
warning(sprintf(c('%i parameters did not converge in numerical SEM derivative.\n',
'Try using different starting values or passing GenRandomPars=TRUE'),
sum(!SEMconv)),
call.=FALSE)
SEMconv <- FALSE
} else SEMconv <- TRUE
DM <- do.call(rbind, DM)
if(!opts$technical$parallel)
.mirtClusterEnv$ncores <- ncores
ESTIMATE$pars <- reloadPars(longpars=ESTIMATE$longpars, pars=ESTIMATE$pars,
ngroups=Data$ngroups, J=Data$nitems)
DM[, is.latent] <- DM[is.latent, ]
DM[is.latent, is.latent] <- 0
info <- try(solve(-solve(ESTIMATE$hess) %*% solve(diag(ncol(DM)) - DM)), silent=TRUE)
info[,is.latent] <- t(info[is.latent, ,drop=FALSE])
if(opts$technical$symmetric) info <- (info + t(info)) / 2
if(is(info, 'try-error')){
if(opts$warn)
warning('Information matrix in SEM could not be computed due to instability',
call.=FALSE)
} else ESTIMATE <- loadESTIMATEinfo(info=info, ESTIMATE=ESTIMATE, constrain=constrain,
warn=opts$warn)
}
} else if(opts$SE.type == 'numerical' && opts$method == 'BL'){
ESTIMATE <- loadESTIMATEinfo(info=-ESTIMATE$hess, ESTIMATE=ESTIMATE, constrain=constrain,
warn=opts$warn)
} else if(opts$SE.type %in% c('Richardson', 'forward', 'central') &&
!(opts$method %in% c('MHRM', 'SEM', 'MIXED'))){
ESTIMATE <- SE.Numerical(pars=ESTIMATE$pars, Theta=ESTIMATE$Theta, theta=theta, PrepList=PrepList, Data=Data,
dentype=opts$dentype, itemloc=PrepList[[1L]]$itemloc, ESTIMATE=ESTIMATE,
constrain=constrain, Ls=Ls, specific=oldmodel, sitems=sitems,
CUSTOM.IND=CUSTOM.IND, EHPrior=ESTIMATE$Prior, warn=opts$warn, type=opts$SE.type,
delta=opts$delta, lrPars=ESTIMATE$lrPars, omp_threads=opts$omp_threads)
} else if(opts$SE.type %in% c('MHRM', 'FMHRM') && opts$method == 'EM'){
if(opts$dentype %in% c('EH', 'EHW'))
stop('MHRM standard error methods not available when using empirical histograms', call.=FALSE)
ESTIMATE <- MHRM.group(pars=pars, constrain=constrain, Ls=Ls, PrepList=PrepList, Data=Data,
list = list(NCYCLES=1000L, BURNIN=1L, SEMCYCLES=opts$SEMCYCLES,
KDRAWS=opts$KDRAWS, MHDRAWS=opts$MHDRAWS,
TOL=opts$SEtol, SE=TRUE, SE.type=opts$SE.type,
gain=opts$gain, nfactNames=PrepList[[1L]]$nfactNames,
itemloc=PrepList[[1L]]$itemloc,
nfact=nfact, constrain=constrain, verbose=FALSE, expl=FALSE,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND, message=opts$message,
startlongpars=startlongpars, SE=opts$SE, warn=opts$warn,
plausible.draws=0L, MSTEPTOL=opts$MSTEPTOL, Moptim='NR1',
keep_vcov_PD=opts$keep_vcov_PD),
DERIV=DERIV, solnp_args=opts$solnp_args, control=control)
} else if(any(opts$SE.type %in% c('crossprod', 'Louis', 'sandwich.Louis', 'sandwich')) &&
!(opts$method %in% c('MHRM', 'SEM', 'MIXED'))){
if(opts$method %in% c('QMCEM', 'MCEM')){
if(opts$warn)
warning(sprintf('SE.type not supported when using method = \"%s\"', opts$method),
call. = FALSE)
} else {
ESTIMATE <- SE.simple(PrepList=PrepList, ESTIMATE=ESTIMATE, Theta=Theta, Data=Data,
constrain=constrain, Ls=Ls, N=nrow(data), type=opts$SE.type,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND, warn=opts$warn,
message=opts$message, complete=ESTIMATE$hess)
}
} else if(opts$SE.type == 'Fisher' && !(opts$method %in% c('MHRM', 'SEM', 'MIXED'))){
if(logPrior != 0 && opts$warn)
warning('Information matrix with the Fisher method does not
account for prior parameter distribution information')
ESTIMATE <- SE.Fisher(PrepList=PrepList, ESTIMATE=ESTIMATE, Theta=Theta, Data=Data,
constrain=constrain, Ls=Ls, full=opts$full,
CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND, warn=opts$warn,
omp_threads=opts$omp_threads)
}
ESTIMATE$cycles <- tmp$cycles
ESTIMATE$Prior <- tmp$Prior
ESTIMATE$Etable <- tmp$Etable
rm(tmp)
}
opts$times$end.time.SE <- proc.time()[3L]
opts$times$start.time.post <- proc.time()[3L]
cmods <- vector('list', Data$ngroups)
if(is.null(opts$theta_lim))
opts$theta_lim <- numeric(1)
lrPars <- ESTIMATE$lrPars
class(lrPars) <- 'S4'
for(g in seq_len(Data$ngroups)){
if(opts$method == 'MIXED' || opts$dentype == "discrete"){
F <- matrix(NA)
h2 <- numeric(1)
} else {
F <- Lambdas(ESTIMATE$pars[[g]], Names=colnames(data))
colnames(F) <- PrepList[[1L]]$factorNames
h2 <- rowSums(F^2)
}
cmods[[g]] <- new('SingleGroupClass', ParObjects=list(pars=ESTIMATE$pars[[g]], lrPars=lrPars,
random=ESTIMATE$random,
lr.random=ESTIMATE$lr.random),
Data = list(K=Data$K, nitems=Data$nitems),
Model = list(itemloc=PrepList[[1L]]$itemloc, nfact=nfact, constrain=constrain,
factorNames=PrepList[[1L]]$factorNames,
itemtype=PrepList[[1L]]$itemtype,
prodlist=PrepList[[1L]]$prodlist),
Options = list(method = 'MHRM', exploratory=PrepList[[1L]]$exploratory,
theta_lim=opts$theta_lim, dentype=opts$dentype),
Fit = list(G2=G2group[g], F=F, h2=h2),
Internals = list(Pl = rlist[[g]]$expected, CUSTOM.IND=CUSTOM.IND,
SLOW.IND=SLOW.IND))
if(opts$dentype %in% c("discrete", 'EH', 'EHW', 'Davidian', 'custom')){
cmods[[g]]@Model$Theta <- Theta
cmods[[g]]@Internals$Prior <- list(ESTIMATE$Prior[[g]])
}
}
#missing stats for MHRM
if(opts$method %in% c('MHRM', 'MIXED', 'SEM') &&
(!opts$logLik_if_converged || !(!ESTIMATE$converge && opts$logLik_if_converged))){
logLik <- G2 <- SElogLik <- 0
if(opts$draws > 0L){
if(opts$verbose) cat("\nCalculating log-likelihood...\n")
flush.console()
if(!opts$technical$parallel){
ncores <- .mirtClusterEnv$ncores
.mirtClusterEnv$ncores <- 1L
}
for(g in seq_len(Data$ngroups)){
cmods[[g]]@Data <- list(data=Data$data[Data$group == Data$groupName[g], ],
fulldata=Data$fulldata[[g]], tabdata=Data$tabdata,
Freq=list(Data$Freq[[g]]), K=Data$K)
cmods[[g]] <- calcLogLik(cmods[[g]], opts$draws, G2 = 'return',
lrPars=ESTIMATE$lrPars)
cmods[[g]]@Data <- list(K=Data$K, nitems=Data$nitems)
logLik <- logLik + cmods[[g]]@Fit$logLik
logPrior <- logPrior + cmods[[g]]@Fit$logPrior
SElogLik <- SElogLik + cmods[[g]]@Fit$SElogLik
G2 <- G2 + cmods[[g]]@Fit$G2
Pl[[g]] <- cmods[[g]]@Internals$Pl
}
if(!opts$technical$parallel)
.mirtClusterEnv$ncores <- ncores
}
}
####post estimation stats
if(opts$Moptim %in% c('solnp', 'nloptr')){
if(!is.null(opts$solnp_args$eqfun))
df <- df + length(opts$solnp_args$eqfun(ESTIMATE$shortpars, list()))
if(!is.null(opts$solnp_args$eval_g_eq))
df <- df + length(opts$solnp_args$eval_g_eq(ESTIMATE$shortpars, list()))
}
r <- rr
N <- sum(r)
tmp <- dfsubtr
AIC <- (-2) * logLik + 2 * tmp
BIC <- (-2) * logLik + tmp*log(N)
SABIC <- (-2) * logLik + tmp*log((N+2)/24)
HQ <- (-2) * logLik + 2*tmp*log(log(N))
p.G2 <- 1 - pchisq(G2,df)
RMSEA.G2 <- rmsea(X2=G2, df=df, N=N)
null.mod <- unclass(new('SingleGroupClass'))
TLI.G2 <- CFI.G2 <- NaN
if(opts$calcNull && length(r) * 3L < prod(Data$K) && opts$warn)
warning(c('Full table of responses is very sparse. ',
'Goodness-of-fit statistics may be very inaccurate'), call.=FALSE)
if(!opts$NULL.MODEL && opts$method != 'MIXED' && opts$calcNull && nmissingtabdata == 0L){
null.mod <- try(unclass(computeNullModel(data=data, key=key,
group=if(length(pars) > 1L) group else NULL)))
if(is(null.mod, 'try-error')){
if(opts$warn)
warning('Null model calculation did not converge.')
null.mod <- unclass(new('SingleGroupClass'))
} else if(!is.nan(G2)) {
TLI.G2 <- tli(X2=G2, X2.null=null.mod@Fit$G2, df=df, df.null=null.mod@Fit$df)
CFI.G2 <- cfi(X2=G2, X2.null=null.mod@Fit$G2, df=df, df.null=null.mod@Fit$df)
}
}
if(nmissingtabdata > 0L)
p.G2 <- RMSEA.G2 <- G2 <- TLI.G2 <- CFI.G2 <- NaN
if(is.null(parprior)) parprior <- list()
if(is.null(opts$quadpts)) opts$quadpts <- NaN
opts$times$end.time.post <- proc.time()[3L]
time <- opts$times
opts$times <- NULL
# set macro objects
Options <- opts
Options$exploratory <- PrepList[[1L]]$exploratory
Fit <- list(G2=G2, p=p.G2, TLI=TLI.G2, CFI=CFI.G2, RMSEA=RMSEA.G2, df=df,
AIC=AIC, BIC=BIC, SABIC=SABIC, HQ=HQ, logLik=logLik,
logPrior=logPrior, SElogLik=SElogLik, F=F, h2=h2)
pis <- if(opts$dentype == 'mixture')
ExtractMixtures(lapply(cmods, function(x) x@ParObjects$pars)) else NULL
Model <- list(model=oldmodel, factorNames=PrepList[[1L]]$factorNames, itemtype=PrepList[[1L]]$itemtype,
itemloc=PrepList[[1L]]$itemloc, nfact=nfact, pis=pis,
Theta=Theta, constrain=constrain, parprior=parprior, nest=as.integer(dfsubtr),
invariance=invariance, lrPars=lrPars, formulas=attr(mixed.design, 'formula'),
prodlist=PrepList[[1L]]$prodlist, nestpars=nestpars)
if(!is.null(opts$technical$Etable)){
Model$Etable <- ESTIMATE$rlist
Model$Etable$Theta <- Theta
}
Data$covdata <- if(length(lrPars)) lrPars@df else attr(mixed.design, 'covdata')
Data$itemdesign <- attr(mixed.design, 'itemdesign')
ParObjects <- list(pars=cmods, lrPars=lrPars, random=ESTIMATE$random, lr.random=ESTIMATE$lr.random)
OptimInfo <- list(iter=ESTIMATE$cycles, converged=ESTIMATE$converge, cand.t.var=ESTIMATE$cand.t.var,
condnum=NA, secondordertest=NA, SEMconv=SEMconv, aveAR=ESTIMATE$aveAR)
vcov <- matrix(NA, 1, 1)
if(Options$SE){
information <- ESTIMATE$info
if(!is.null(opts$technical$infoAsVcov)){
vcov <- information
} else {
if(!ESTIMATE$fail_invert_info){
isna <- is.na(diag(information))
info <- information[!isna, !isna]
vcov <- matrix(NA, ncol(information), ncol(information))
rownames(vcov) <- colnames(vcov) <- colnames(information)
vcov2 <- try(solve(info), silent=TRUE)
vcov[!isna, !isna] <- vcov2
if(!is(vcov2, 'try-error')){
OptimInfo$condnum <- kappa(info, exact=TRUE)
OptimInfo$secondordertest <- secondOrderTest(info)
} else OptimInfo$secondordertest <- FALSE
} else {
OptimInfo$secondordertest <- FALSE
}
}
} else OptimInfo$secondordertest <- NA
Internals <- list(collectLL=ESTIMATE$collectLL, Prior=ESTIMATE$Prior, Pl=Pl,
shortpars=as.numeric(ESTIMATE$shortpars), key=key,
bfactor=list(), CUSTOM.IND=CUSTOM.IND, SLOW.IND=SLOW.IND,
survey.weights=survey.weights, theta_lim = opts$theta_lim,
customGroup=customGroup, customItems=customItems)
if(opts$method == 'EM'){
tmp <- lapply(ESTIMATE$Etable, function(tab)
data.frame(Theta, posterior=rowSums(tab$r1)))
if(length(tmp)){
names(tmp) <- Data$groupNames
Internals$thetaPosterior <- tmp
}
}
if(opts$storeEtable)
Internals$Etable <- ESTIMATE$Etable
if(opts$storeEMhistory)
Internals$EMhistory <- ESTIMATE$EMhistory
if(opts$method == 'SEM') Options$TOL <- NA
if(opts$odentype == "discrete"){
Fit$F <- Fit$h2 <- NULL
mod <- new('DiscreteClass',
Data=Data,
Options=Options,
Fit=Fit,
Model=Model,
ParObjects=ParObjects,
OptimInfo=OptimInfo,
Internals=Internals,
vcov=vcov)
} else {
if(Data$ngroups == 1L){
ParObjects$pars <- cmods[[1L]]@ParObjects$pars
Internals$Pl <- Internals$Pl[[1L]]
if(opts$method == 'MIXED'){
Fit$p <- NaN
mod <- new('MixedClass',
Data=Data,
Options=Options,
Fit=Fit,
Model=Model,
ParObjects=ParObjects,
OptimInfo=OptimInfo,
Internals=Internals,
vcov=vcov)
} else {
if(Options$exploratory){
FF <- F %*% t(F)
V <- eigen(FF)$vector[ ,1L:nfact]
L <- eigen(FF)$values[1L:nfact]
if (nfact == 1L) F <- as.matrix(V * sqrt(L))
else F <- V %*% sqrt(diag(L))
if (sum(F[ ,1L] < 0)) F <- (-1) * F
colnames(F) <- paste("F", 1L:ncol(F), sep="")
h2 <- rowSums(F^2)
} else {
if(opts$method == 'EM')
Internals$bfactor <- list(prior=ESTIMATE$prior,
Priorbetween=ESTIMATE$Priorbetween,
sitems=ESTIMATE$sitems, specific=specific)
}
mod <- new('SingleGroupClass',
Data=Data,
Options=Options,
Fit=Fit,
Model=Model,
ParObjects=ParObjects,
OptimInfo=OptimInfo,
Internals=Internals,
vcov=vcov)
}
} else {
if(opts$method == 'EM')
Internals$bfactor <- list(prior=ESTIMATE$prior,
Priorbetween=ESTIMATE$Priorbetween,
sitems=ESTIMATE$sitems, specific=specific)
cls <- ifelse(opts$odentype == 'mixture', 'MixtureClass', 'MultipleGroupClass')
mod <- new(cls,
Data=Data,
Options=Options,
Fit=Fit,
Model=Model,
ParObjects=ParObjects,
OptimInfo=OptimInfo,
Internals=Internals,
vcov=vcov)
}
}
mod@time <- c("TOTAL:" = as.numeric(proc.time()[3L] - time$start.time),
Data = as.numeric(time$end.time.Data - time$start.time.Data),
ESTIMATE$time,
SE = as.numeric(time$end.time.SE - time$start.time.SE),
Post = as.numeric(time$end.time.post - time$start.time.post))
return(mod)
}
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.