# theta combinations
thetaComb <- function(theta, nfact)
{
if (nfact == 1L){
Theta <- matrix(theta)
} else {
thetalist <- vector('list', nfact)
for(i in seq_len(nfact))
thetalist[[i]] <- theta
Theta <- as.matrix(expand.grid(thetalist))
}
return(Theta)
}
# Product terms
prodterms <- function(theta0, prodlist)
{
products <- matrix(1, ncol = length(prodlist), nrow = nrow(theta0))
for(i in seq_len(length(prodlist))){
tmp <- prodlist[[i]]
for(j in 1L:length(tmp))
products[ ,i] <- products[ ,i] * theta0[ ,tmp[j]]
}
ret <- cbind(theta0,products)
ret
}
# MH sampler for theta values
draw.thetas <- function(theta0, pars, fulldata, itemloc, cand.t.var, prior.t.var,
prior.mu, prodlist, OffTerm, CUSTOM.IND)
{
makedraws <- try({
N <- nrow(fulldata)
unif <- runif(N)
sigma <- if(ncol(theta0) == 1L) matrix(cand.t.var) else diag(rep(cand.t.var,ncol(theta0)))
total_0 <- attr(theta0, 'log.lik_full')
theta1prod <- theta1 <- theta0 + mirt_rmvnorm(N, sigma = sigma)
if(is.null(total_0)) theta1prod <- theta1 <- theta0 #for intial draw
if(length(prodlist))
theta1prod <- prodterms(theta1, prodlist)
total_1 <- complete.LL(theta=theta1prod, pars=pars, nfact=ncol(theta1), prior.mu=prior.mu,
prior.t.var=prior.t.var, OffTerm=OffTerm,
CUSTOM.IND=CUSTOM.IND, itemloc=itemloc, fulldata=fulldata)
if(is.null(total_0)){ #for intial draw
attr(theta1, 'log.lik_full') <- total_1
return(theta1)
}
diff <- total_1 - total_0
accept <- unif < exp(diff)
theta1[!accept, ] <- theta0[!accept, ]
total_1[!accept] <- total_0[!accept]
log.lik <- sum(total_1)
attr(theta1, "Proportion Accepted") <- sum(accept)/N
attr(theta1, "log.lik") <- log.lik
attr(theta1, 'log.lik_full') <- total_1
}, silent = TRUE)
if(is(makedraws, 'try-error'))
stop('MH sampler failed. Model is likely unstable or may need better starting values',
.call=FALSE)
return(theta1)
}
complete.LL <- function(theta, pars, nfact, prior.mu, prior.t.var,
OffTerm, CUSTOM.IND, itemloc, fulldata){
log_den <- mirt_dmvnorm(theta[,seq_len(nfact), drop=FALSE], prior.mu, prior.t.var, log=TRUE)
itemtrace <- computeItemtrace(pars=pars, Theta=theta, itemloc=itemloc,
offterm=OffTerm, CUSTOM.IND=CUSTOM.IND)
rowSums(fulldata * log(itemtrace)) + log_den
}
imputePars <- function(pars, imputenums, constrain, pre.ev){
shift <- mirt_rmvnorm(1L, mean=numeric(length(pre.ev$values)), pre.ev=pre.ev)
for(i in seq_len(length(pars))){
pn <- pars[[i]]@parnum
pick2 <- imputenums %in% pn
pick1 <- pn %in% imputenums
pars[[i]]@par[pick1] <- pars[[i]]@par[pick1] + shift[pick2]
if(is(pars[[i]], 'graded')){
where <- (length(pars[[i]]@par) - pars[[i]]@ncat + 2L):length(pars[[i]]@par)
if(!(all(sort(pars[[i]]@par[where], decreasing=TRUE) == pars[[i]]@par[where])))
stop('Drawn values out of order', call.=FALSE)
} else if(is(pars[[i]], 'grsm')){
where <- (length(pars[[i]]@par) - pars[[i]]@ncat + 1L):(length(pars[[i]]@par)-1L)
if(!(all(sort(pars[[i]]@par[where], decreasing=TRUE) == pars[[i]]@par[where])))
stop('Drawn values out of order', call.=FALSE)
}
}
for(con in seq_len(length(constrain))){
tmp <- shift[imputenums %in% constrain[[con]][1L]]
if(length(tmp)){
for(i in seq_len(length(pars))){
pick <- pars[[i]]@parnum %in% constrain[[con]][-1L]
pars[[i]]@par[pick] <- tmp + pars[[i]]@par[pick]
}
}
}
return(pars)
}
imputePars2 <- function(MGmod, shortpars, longpars, imputenums, pre.ev){
while(TRUE){
shift <- mirt_rmvnorm(1L, mean=shortpars, pre.ev=pre.ev)
longpars[imputenums] <- shift[1L,]
constrain <- MGmod@Model$constrain
longpars <- longpars_constrain(longpars=longpars, constrain=constrain)
pars <- list(MGmod@ParObjects$pars[[1L]]@ParObjects$pars, MGmod@ParObjects$pars[[2L]]@ParObjects$pars)
pars <- reloadPars(longpars=longpars, pars=pars, ngroups=2L, J=length(pars[[1L]])-1L)
if(any(MGmod@Model$itemtype %in% c('graded', 'grsm'))){
pick <- c(MGmod@Model$itemtype %in% c('graded', 'grsm'), FALSE)
if(!all(sapply(pars[[1L]][pick], CheckIntercepts) &
sapply(pars[[2L]][pick], CheckIntercepts))) next
}
break
}
MGmod@ParObjects$pars[[1L]]@ParObjects$pars <- pars[[1L]]
MGmod@ParObjects$pars[[2L]]@ParObjects$pars <- pars[[2L]]
MGmod
}
# Rotation function
Rotate <- function(F, rotate, Target = NULL, par.strip.text = NULL, par.settings = NULL, ...)
{
if(ncol(F) == 1L) rotF <- list()
if(rotate == 'none') rotF <- list(loadings=F, Phi=diag(ncol(F)), orthogonal=TRUE)
if(rotate == 'promax'){
mypromax <- function (x, m = 4) {
#borrowed and modified from stats::promax on Febuary 13, 2013
if (ncol(x) < 2L)
return(x)
dn <- dimnames(x)
xx <- varimax(x)
x <- xx$loadings
Q <- x * abs(x)^(m - 1)
U <- lm.fit(x, Q)$coefficients
d <- diag(solve(t(U) %*% U))
U <- U %*% diag(sqrt(d))
dimnames(U) <- NULL
z <- x %*% U
U <- xx$rotmat %*% U
ui <- solve(U)
Phi <- ui %*% t(ui)
dimnames(z) <- dn
class(z) <- "loadings"
list(loadings = z, rotmat = U, Phi = Phi, orthogonal = FALSE)
}
rotF <- mypromax(F, ...)
}
if(rotate == 'oblimin') rotF <- GPArotation::oblimin(F, ...)
if(rotate == 'quartimin') rotF <- GPArotation::quartimin(F, ...)
if(rotate == 'targetT') rotF <- GPArotation::targetT(F, Target = Target, ...)
if(rotate == 'targetQ') rotF <- GPArotation::targetQ(F, Target = Target, ...)
if(rotate == 'pstT') rotF <- GPArotation::pstT(F, Target = Target, ...)
if(rotate == 'pstQ') rotF <- GPArotation::pstQ(F, Target = Target, ...)
if(rotate == 'oblimax') rotF <- GPArotation::oblimax(F, ...)
if(rotate == 'entropy') rotF <- GPArotation::entropy(F, ...)
if(rotate == 'quartimax') rotF <- GPArotation::quartimax(F, ...)
if(rotate == 'varimax') rotF <- GPArotation::Varimax(F, ...)
if(rotate == 'simplimax') rotF <- GPArotation::simplimax(F, ...)
if(rotate == 'bentlerT') rotF <- GPArotation::bentlerT(F, ...)
if(rotate == 'bentlerQ') rotF <- GPArotation::bentlerQ(F, ...)
if(rotate == 'tandemI') rotF <- GPArotation::tandemI(F, ...)
if(rotate == 'tandemII') rotF <- GPArotation::tandemII(F, ...)
if(rotate == 'geominT') rotF <- GPArotation::geominT(F, ...)
if(rotate == 'geominQ') rotF <- GPArotation::geominQ(F, ...)
if(rotate == 'cfT') rotF <- GPArotation::cfT(F, ...)
if(rotate == 'cfQ') rotF <- GPArotation::cfQ(F, ...)
if(rotate == 'infomaxT') rotF <- GPArotation::infomaxT(F, ...)
if(rotate == 'infomaxQ') rotF <- GPArotation::infomaxQ(F, ...)
if(rotate == 'mccammon') rotF <- GPArotation::mccammon(F, ...)
if(rotate == 'bifactorT') rotF <- GPArotation::bifactorT(F, ...)
if(rotate == 'bifactorQ') rotF <- GPArotation::bifactorQ(F, ...)
return(unclass(rotF))
}
# Gamma correlation, mainly for obtaining a sign
gamma.cor <- function(x)
{
concordant <- function(x){
mat.lr <- function(r, c, r.x, c.x){
lr <- x[(r.x > r) & (c.x > c)]
sum(lr)
}
r.x <- row(x)
c.x <- col(x)
sum(x * mapply(mat.lr, r = r.x, c = c.x, MoreArgs=list(r.x=r.x, c.x=c.x)))
}
discordant <- function(x){
mat.ll <- function(r, c, r.x, c.x){
ll <- x[(r.x > r) & (c.x < c)]
sum(ll)
}
r.x <- row(x)
c.x <- col(x)
sum(x * mapply(mat.ll, r = r.x, c = c.x, MoreArgs=list(r.x=r.x, c.x=c.x)))
}
c <- concordant(x)
d <- discordant(x)
gamma <- (c - d) / (c + d)
gamma
}
# Approximation to polychoric matrix for initial values
cormod <- function(fulldata, K, guess, smooth = TRUE, use = 'pairwise.complete.obs')
{
fulldata <- as.matrix(fulldata)
cormat <- suppressWarnings(cor(fulldata, use=use))
diag(cormat) <- 1
cormat[is.na(cormat)] <- 0
cormat <- abs(cormat)^(1/1.15) * sign(cormat)
if(smooth)
cormat <- smooth.cor(cormat)
cormat
}
# Rotate lambda coefficients
rotateLambdas <- function(so){
F <- so$rotF
h2 <- so$h2
h <- matrix(rep(sqrt(1 - h2), ncol(F)), ncol = ncol(F))
a <- F / h
a
}
d2r <-function(d) pi*d/180
closeEnough <- function(x, low, up) all(x >= low & x <= up)
logit <- function(x){
ret <- qlogis(x)
ret <- ifelse(x == 0, -999, ret)
ret <- ifelse(x == 1, 999, ret)
ret
}
antilogit <- function(x) plogis(x)
MPinv <- function(mat){
svd <- svd(mat)
d <- svd$d; v <- svd$v; u <- svd$u
P <- d > max(sqrt(.Machine$double.eps) * d[1L], 0)
if(all(P)){
mat <- v %*% (1/d * t(u))
} else {
mat <- v[ , P, drop=FALSE] %*% ((1/d[P]) * t(u[ , P, drop=FALSE]))
}
mat
}
test_info <- function(pars, Theta, Alist, K){
infolist <- list()
for(cut in seq_len(length(Alist))){
A <- Alist[[cut]]
info <- rep(0,nrow(Theta))
for(j in seq_len(length(K))){
info <- info + ItemInfo(pars[[j]], A[j,], Theta)
}
infolist[[cut]] <- info
}
tmp <- 0
for(i in seq_len(length(infolist))){
tmp <- tmp + infolist[[i]]
}
info <- tmp/length(infolist)
info
}
Lambdas <- function(pars, Names){
J <- length(pars) - 1L
lambdas <- matrix(NA, J, length(ExtractLambdas(pars[[1L]])))
gcov <- ExtractGroupPars(pars[[J+1L]])$gcov
if(ncol(gcov) < ncol(lambdas)){
tmpcov <- diag(ncol(lambdas))
tmpcov[seq_len(ncol(gcov)), seq_len(ncol(gcov))] <- gcov
gcov <- tmpcov
}
rownames(lambdas) <- Names
for(i in seq_len(J)){
tmp <- pars[[i]]
lambdas[i,] <- ExtractLambdas(tmp) /1.702
}
norm <- sqrt(1 + rowSums(lambdas^2))
dcov <- if(ncol(gcov) > 1L) diag(sqrt(diag(gcov))) else matrix(sqrt(diag(gcov)))
F <- as.matrix(lambdas/norm) %*% dcov
F
}
#change long pars for groups into mean in sigma
ExtractGroupPars <- function(x){
if(x@itemclass < 0L) return(list(gmeans=0, gcov=matrix(1)))
nfact <- x@nfact
gmeans <- x@par[seq_len(nfact)]
tmp <- x@par[-seq_len(nfact)]
gcov <- matrix(0, nfact, nfact)
gcov[lower.tri(gcov, diag=TRUE)] <- tmp
if(nfact != 1L)
gcov <- gcov + t(gcov) - diag(diag(gcov))
return(list(gmeans=gmeans, gcov=gcov))
}
reloadConstr <- function(par, constr, obj){
par2 <- rep(NA, length(constr[[1L]]))
notconstr <- rep(TRUE, length(par2))
for(i in seq_len(length(constr))){
par2[constr[[i]]] <- par[i]
notconstr[constr[[i]]] <- FALSE
}
par2[notconstr] <- par[(length(constr)+1L):length(par)]
ind <- 1L
for(i in seq_len(length(obj))){
obj[[i]]@par[obj[[i]]@est] <- par2[ind:(ind + sum(obj[[i]]@est) - 1L)]
ind <- ind + sum(obj[[i]]@est)
}
return(obj)
}
bfactor2mod <- function(model, J){
tmp <- tempfile('tempfile')
unique <- sort(unique(model))
index <- seq_len(J)
tmp2 <- c()
for(i in seq_len(length(unique))){
ind <- na.omit(index[model == unique[i]])
comma <- rep(',', 2*length(ind))
TF <- rep(c(TRUE,FALSE), length(ind))
comma[TF] <- ind
comma[length(comma)] <- ""
tmp2 <- c(tmp2, c(paste('\nS', i, ' =', sep=''), comma))
}
cat(tmp2, file=tmp)
model <- mirt.model(file=tmp, quiet = TRUE)
unlink(tmp)
return(model)
}
updateTheta <- function(npts, nfact, pars, QMC = FALSE){
ngroups <- length(pars)
pick <- length(pars[[1L]])
Theta <- vector('list', ngroups)
for(g in seq_len(ngroups)){
gp <- ExtractGroupPars(pars[[g]][[pick]])
theta <- if(QMC){
QMC_quad(npts=npts, nfact=nfact, lim = c(0,1))
} else {
MC_quad(npts=npts, nfact=nfact, lim = c(0,1))
}
Theta[[g]] <- t(t(theta) + gp$gmeans) %*% t(chol(gp$gcov))
}
Theta
}
updatePrior <- function(pars, gTheta, list, ngroups, nfact, J,
dentype, sitems, cycles, rlist, lrPars = list(), full=FALSE,
MC = FALSE){
prior <- Prior <- Priorbetween <- vector('list', ngroups)
if(dentype == 'EH'){
Prior[[1L]] <- list$EHPrior[[1L]]
} else if(dentype == 'custom'){
for(g in seq_len(ngroups)){
gp <- pars[[g]][[J+1L]]
Prior[[g]] <- gp@den(gp, gTheta[[g]])
Prior[[g]] <- Prior[[g]] / sum(Prior[[g]])
}
} else if(dentype == 'discrete'){
for(g in seq_len(ngroups)){
gp <- pars[[g]][[J+1L]]
if(full){
Prior[[g]] <- gp@den(gp, gTheta[[g]], mus=lrPars@mus)
Prior[[g]] <- Prior[[g]]/rowSums(Prior[[g]])
} else {
Prior[[g]] <- gp@den(gp, gTheta[[g]])
Prior[[g]] <- Prior[[g]] / sum(Prior[[g]])
}
}
} else {
for(g in seq_len(ngroups)){
gp <- ExtractGroupPars(pars[[g]][[J+1L]])
if(dentype == 'bfactor'){
theta <- pars[[g]][[J+1L]]@theta
Thetabetween <- pars[[g]][[J+1L]]@Thetabetween
p <- matrix(0, nrow(gTheta[[g]]), ncol(sitems))
pp <- matrix(0, nrow(theta), ncol(sitems))
for(i in seq_len(ncol(sitems))){
sel <- c(seq_len(nfact-ncol(sitems)), i + nfact - ncol(sitems))
p[,i] <- mirt_dmvnorm(gTheta[[g]][ ,sel], gp$gmeans[sel], gp$gcov[sel,sel,drop=FALSE])
pp[,i] <- dnorm(theta, gp$gmeans[sel[length(sel)]],
sqrt(gp$gcov[sel[length(sel)],sel[length(sel)],drop=FALSE]))
}
pb <- mirt_dmvnorm(Thetabetween, gp$gmeans[seq_len(ncol(Thetabetween))],
gp$gcov[seq_len(ncol(Thetabetween)), seq_len(ncol(Thetabetween)), drop=FALSE])
Priorbetween[[g]] <- pb / sum(pb)
Prior[[g]] <- t(t(p) / colSums(p))
prior[[g]] <- t(t(pp) / colSums(pp))
next
}
if(full){
Prior[[g]] <- mirt_dmvnorm(gTheta[[g]][ ,seq_len(nfact),drop=FALSE], lrPars@mus, gp$gcov,
quad=TRUE)
Prior[[g]] <- Prior[[g]]/rowSums(Prior[[g]])
} else {
Prior[[g]] <- mirt_dmvnorm(gTheta[[g]][ ,seq_len(nfact),drop=FALSE], gp$gmeans, gp$gcov)
Prior[[g]] <- Prior[[g]]/sum(Prior[[g]])
}
}
}
if(dentype == 'EH'){
if(cycles > 1L){
for(g in seq_len(ngroups))
Prior[[g]] <- rowSums(rlist[[g]][[1L]]) / sum(rlist[[g]][[1L]])
} else {
for(g in seq_len(ngroups)){
Prior[[g]] <- mirt_dmvnorm(gTheta[[g]], 0, matrix(1))
Prior[[g]] <- Prior[[g]]/sum(Prior[[g]])
}
}
} else if(!is.null(list$customPriorFun)){
for(g in seq_len(ngroups)){
Prior[[g]] <- list$customPriorFun(gTheta[[g]], Etable=rlist[[g]][[1L]])
Prior[[g]] <- Prior[[g]]/sum(Prior[[g]])
}
}
if(MC){
if(full){
for(g in seq_len(ngroups))
Prior[[g]] <- matrix(rep(1 / length(gTheta[[g]])),
nrow(lrPars@mus), nrow(gTheta[[g]]))
} else {
for(g in seq_len(ngroups))
Prior[[g]] <- matrix(rep(1 / length(Prior[[g]]), length(Prior[[g]])))
}
}
return(list(prior=prior, Prior=Prior, Priorbetween=Priorbetween))
}
UpdateConstrain <- function(pars, constrain, invariance, nfact, nLambdas, J, ngroups, PrepList,
method, itemnames, model, groupNames)
{
if(!is.numeric(model[[1L]])){
if(any(model[[1L]]$x[,1L] == 'CONSTRAIN')){
groupNames <- as.character(groupNames)
names(pars) <- groupNames
input <- model[[1L]]$x[model[[1L]]$x[,1L] == 'CONSTRAIN', 2L]
input <- gsub(' ', replacement='', x=input)
elements <- strsplit(input, '\\),\\(')[[1L]]
elements <- gsub('\\(', replacement='', x=elements)
elements <- gsub('\\)', replacement='', x=elements)
esplit <- strsplit(elements, ',')
esplit <- lapply(esplit, function(x, groupNames)
if(!(x[length(x)] %in% c(groupNames, 'all'))) c(x, 'all') else x,
groupNames=as.character(groupNames))
esplit <- lapply(esplit, function(x){
newx <- c()
if(length(x) < 3L)
stop('CONTRAIN = ... has not been supplied enough arguments', call.=FALSE)
for(i in seq_len(length(x)-2L)){
if(grepl('-', x[i])){
tmp <- as.numeric(strsplit(x[i], '-')[[1L]])
newx <- c(newx, tmp[1L]:tmp[2L])
} else newx <- c(newx, x[i])
}
x <- c(newx, x[length(x)-1L], x[length(x)])
x
})
for(i in seq_len(length(esplit))){
if(!(esplit[[i]][length(esplit[[i]])] %in% c(groupNames, 'all')))
stop('Invalid group name passed to CONSTRAIN = ... syntax.', call.=FALSE)
if(esplit[[i]][length(esplit[[i]])] == 'all'){
for(g in seq_len(ngroups)){
constr <- c()
p <- pars[[g]]
sel <- suppressWarnings(
as.numeric(esplit[[i]][seq_len(length(esplit[[i]])-1L)]))
picknames <- c(is.na(sel), FALSE)
sel <- na.omit(sel)
if(length(sel) == 1L){
for(j in seq_len(length(sel))){
pick <- p[[sel[j]]]@parnum[names(p[[sel[j]]]@est) %in%
esplit[[i]][picknames]]
if(!length(pick))
stop('CONSTRAIN = ... indexed a parameter that was not relavent for item ', sel[j],
call.=FALSE)
constr <- c(constr, pick)
}
} else {
if(sum(picknames) > 1L){
if(sum(picknames) != length(sel))
stop('Number of items selected not equal to number of parameter names',
call.=FALSE)
constr <- numeric(length(sel))
for(j in seq_len(length(sel))){
whc <- esplit[[i]][which(picknames)[j]]
constr[j] <- p[[sel[j]]]@parnum[names(p[[sel[j]]]@est) == whc]
}
} else {
for(j in seq_len(length(sel))){
pick <- p[[sel[j]]]@parnum[names(p[[sel[j]]]@est) %in%
esplit[[i]][picknames]]
if(!length(pick))
stop('CONSTRAIN = ... indexed a parameter that was not relavent for item ', sel[j],
call.=FALSE)
constr <- c(constr, pick)
}
}
}
constrain[[length(constrain) + 1L]] <- constr
}
} else {
constr <- c()
p <- pars[[esplit[[i]][length(esplit[[i]])]]]
sel <- as.numeric(esplit[[i]][seq_len(length(esplit[[i]])-2L)])
for(j in seq_len(length(sel))){
pick <- p[[sel[j]]]@parnum[names(p[[sel[j]]]@est) ==
esplit[[i]][length(esplit[[i]])-1L]]
if(!length(pick))
stop('CONSTRAIN = ... indexed a parameter that was not relavent for item ', sel[j],
call.=FALSE)
constr <- c(constr, pick)
}
constrain[[length(constrain) + 1L]] <- constr
}
}
}
if(any(model[[1L]]$x[,1L] == 'CONSTRAINB')){
if(length(unique(groupNames)) == 1L)
stop('CONSTRAINB model argument not valid for single group models', call.=FALSE)
groupNames <- as.character(groupNames)
names(pars) <- groupNames
input <- model[[1L]]$x[model[[1L]]$x[,1L] == 'CONSTRAINB', 2L]
input <- gsub(' ', replacement='', x=input)
elements <- strsplit(input, '\\),\\(')[[1L]]
elements <- gsub('\\(', replacement='', x=elements)
elements <- gsub('\\)', replacement='', x=elements)
esplit <- strsplit(elements, ',')
esplit <- lapply(esplit, function(x, groupNames)
if(!(x[length(x)] %in% c(groupNames, 'all'))) c(x, 'all') else x,
groupNames=as.character(groupNames))
esplit <- lapply(esplit, function(x){
newx <- c()
if(length(x) < 3L)
stop('PRIOR = ... has not been supplied enough arguments', call.=FALSE)
for(i in seq_len(length(x)-2L)){
if(grepl('-', x[i])){
tmp <- as.numeric(strsplit(x[i], '-')[[1L]])
newx <- c(newx, tmp[1L]:tmp[2L])
} else newx <- c(newx, x[i])
}
x <- c(newx, x[length(x)-1L], x[length(x)])
x
})
for(i in seq_len(length(esplit))){
sel <- as.numeric(esplit[[i]][seq_len(length(esplit[[i]])-2L)])
for(j in seq_len(length(sel))){
constr <- c()
for(g in seq_len(ngroups)){
p <- pars[[g]]
pick <- p[[sel[j]]]@parnum[names(p[[sel[j]]]@est) ==
esplit[[i]][length(esplit[[i]])-1L]]
if(!length(pick))
stop('CONSTRAINB = ... indexed a parameter that was not relavent accross groups',
call.=FALSE)
constr <- c(constr, pick)
}
constrain[[length(constrain) + 1L]] <- constr
}
}
}
}
#within group item constraints only
for(g in seq_len(ngroups))
for(i in seq_len(length(PrepList[[g]]$constrain)))
constrain[[length(constrain) + 1L]] <- PrepList[[g]]$constrain[[i]]
if('covariances' %in% invariance){ #Fix covariance accross groups (only makes sense with vars = 1)
tmp <- c()
tmpmats <- tmpestmats <- matrix(NA, ngroups, nfact*(nfact+1L)/2)
for(g in seq_len(ngroups)){
tmpmats[g,] <- pars[[g]][[J + 1L]]@parnum[(nfact+1L):length(pars[[g]][[J + 1L]]@parnum)]
tmpestmats[g,] <- pars[[g]][[J + 1L]]@est[(nfact+1L):length(pars[[g]][[J + 1L]]@est)]
}
select <- colSums(tmpestmats) == ngroups
for(i in seq_len(length(select)))
if(select[i])
constrain[[length(constrain) + 1L]] <- tmpmats[seq_len(ngroups), i]
}
if(any(itemnames %in% invariance)){
matched <- na.omit(match(invariance, itemnames))
for(i in matched){
jj <- sum(pars[[1L]][[i]]@est)
stopifnot(jj > 0)
for(j in seq_len(jj)){
tmp <- c()
for(g in seq_len(ngroups))
tmp <- c(tmp, pars[[g]][[i]]@parnum[pars[[g]][[i]]@est][j])
constrain[[length(constrain) + 1L]] <- tmp
}
}
}
if('slopes' %in% invariance){ #Equal factor loadings
tmpmats <- tmpests <- list()
for(g in seq_len(ngroups))
tmpmats[[g]] <- tmpests[[g]] <- matrix(NA, J, nLambdas)
for(g in seq_len(ngroups)){
for(i in seq_len(J)){
tmpmats[[g]][i,] <- pars[[g]][[i]]@parnum[1L:nLambdas]
tmpests[[g]][i,] <- pars[[g]][[i]]@est[1L:nLambdas]
}
}
for(i in seq_len(J)){
for(j in seq_len(nLambdas)){
tmp <- c()
for(g in seq_len(ngroups)){
if(tmpests[[1L]][[i, j]])
tmp <- c(tmp, tmpmats[[g]][i,j])
}
constrain[[length(constrain) + 1L]] <- tmp
}
}
}
if('intercepts' %in% invariance){ #Equal item intercepts (and all other item pars)
tmpmats <- tmpests <- list()
for(g in seq_len(ngroups))
tmpmats[[g]] <- tmpests[[g]] <- list()
for(g in seq_len(ngroups)){
for(i in seq_len(J)){
ind <- (nLambdas+1L):length(pars[[g]][[i]]@parnum)
if(is(pars[[g]][[i]], 'dich')) ind <- ind[seq_len(length(ind)-2L)]
if(is(pars[[g]][[i]], 'partcomp')) ind <- ind[seq_len(length(ind)-1L)]
tmpmats[[g]][[i]] <- pars[[g]][[i]]@parnum[ind]
tmpests[[g]][[i]] <- pars[[g]][[i]]@est[ind]
}
}
for(i in seq_len(J)){
for(j in seq_len(length(tmpmats[[1L]][[i]]))){
tmp <- c()
for(g in seq_len(ngroups)){
if(tmpests[[1L]][[i]][j])
tmp <- c(tmp, tmpmats[[g]][[i]][j])
}
constrain[[length(constrain) + 1L]] <- tmp
}
}
}
#remove redundent constraints
redun <- rep(FALSE, length(constrain))
if(length(constrain)){
for(i in seq_len(length(redun))){
while(TRUE){
lastredun <- redun
for(j in seq_len(length(redun))){
if(i < j && !redun[j] && !redun[i]){
if(any(constrain[[i]] %in% constrain[[j]])){
constrain[[i]] <- unique(c(constrain[[i]], constrain[[j]]))
redun[j] <- TRUE
}
}
}
if(all(lastredun == redun)) break
}
}
}
constrain[redun] <- NULL
return(constrain)
}
expbeta_sv <- function(val1, val2){
ret <- qlogis((val1-1)/(val1 + val2-2))
if(!is.finite(ret)) ret <- qlogis(val1/(val1 + val2))
ret
}
UpdatePrior <- function(PrepList, model, groupNames){
if(!is.numeric(model[[1L]])){
if(!length(model[[1L]]$x[model[[1L]]$x[,1L] == 'PRIOR', 2L])) return(PrepList)
groupNames <- as.character(groupNames)
ngroups <- length(groupNames)
pars <- vector('list', length(PrepList))
for(g in seq_len(length(PrepList)))
pars[[g]] <- PrepList[[g]]$pars
names(pars) <- groupNames
input <- model[[1L]]$x[model[[1L]]$x[,1L] == 'PRIOR', 2L]
input <- gsub(' ', replacement='', x=input)
elements <- strsplit(input, '\\),\\(')[[1L]]
elements <- gsub('\\(', replacement='', x=elements)
elements <- gsub('\\)', replacement='', x=elements)
esplit <- strsplit(elements, ',')
esplit <- lapply(esplit, function(x, groupNames)
if(!(x[length(x)] %in% c(groupNames, 'all'))) c(x, 'all') else x,
groupNames=as.character(groupNames))
esplit <- lapply(esplit, function(x){
newx <- c()
if(length(x) < 5L)
stop('PRIOR = ... has not been supplied enough arguments', call.=FALSE)
for(i in seq_len(length(x)-5L)){
if(grepl('-', x[i])){
tmp <- as.numeric(strsplit(x[i], '-')[[1L]])
newx <- c(newx, tmp[1L]:tmp[2L])
} else newx <- c(newx, x[i])
}
x <- c(newx, x[(length(x)-4L):length(x)])
x
})
for(i in seq_len(length(esplit))){
if(!(esplit[[i]][length(esplit[[i]])] %in% c(groupNames, 'all')))
stop('Invalid group name passed to PRIOR = ... syntax.', call.=FALSE)
if(esplit[[i]][length(esplit[[i]])] == 'all'){
for(g in seq_len(ngroups)){
sel <- as.numeric(esplit[[i]][1L:(length(esplit[[i]])-5L)])
name <- esplit[[i]][length(esplit[[i]])-4L]
type <- esplit[[i]][length(esplit[[i]])-3L]
if(!(type %in% c('norm', 'beta', 'lnorm', 'expbeta')))
stop('Prior type specified in PRIOR = ... not available', call.=FALSE)
type <- switch(type, norm=1L, lnorm=2L, beta=3L, expbeta=4L, 0L)
val1 <- as.numeric(esplit[[i]][length(esplit[[i]])-2L])
val2 <- as.numeric(esplit[[i]][length(esplit[[i]])-1L])
for(j in seq_len(length(sel))){
which <- names(pars[[g]][[sel[j]]]@est) == name
if(!any(which)) stop('Parameter \'', name, '\' does not exist for item ', j,
call.=FALSE)
pars[[g]][[sel[j]]]@any.prior <- TRUE
pars[[g]][[sel[j]]]@prior.type[which] <- type
pars[[g]][[sel[j]]]@prior_1[which] <- val1
pars[[g]][[sel[j]]]@prior_2[which] <- val2
pars[[g]][[sel[j]]]@par[which] <- switch(type,
'1'=val1,
'2'=exp(val1),
'3'=(val1-1)/(val1 + val2 - 2),
'4'=expbeta_sv(val1, val2))
if(type == '2')
pars[[g]][[sel[j]]]@lbound[which] <- 0
if(type == '3'){
pars[[g]][[sel[j]]]@lbound[which] <- 0
pars[[g]][[sel[j]]]@ubound[which] <- 1
}
}
}
} else {
sel <- as.numeric(esplit[[i]][seq_len(length(esplit[[i]])-5L)])
gname <- esplit[[i]][length(esplit[[i]])]
name <- esplit[[i]][length(esplit[[i]])-4L]
type <- esplit[[i]][length(esplit[[i]])-3L]
if(!(type %in% c('norm', 'beta', 'lnorm', 'expbeta')))
stop('Prior type specified in PRIOR = ... not available', call.=FALSE)
type <- switch(type, norm=1L, lnorm=2L, beta=3L, expbeta=4L, 0L)
val1 <- as.numeric(esplit[[i]][length(esplit[[i]])-2L])
val2 <- as.numeric(esplit[[i]][length(esplit[[i]])-1L])
for(j in seq_len(length(sel))){
which <- names(pars[[gname]][[sel[j]]]@est) == name
if(!any(which)) stop('Parameter \'', name, '\' does not exist for item ', j,
call.=FALSE)
pars[[gname]][[sel[j]]]@any.prior <- TRUE
pars[[gname]][[sel[j]]]@prior.type[which] <- type
pars[[gname]][[sel[j]]]@prior_1[which] <- val1
pars[[gname]][[sel[j]]]@prior_2[which] <- val2
}
}
}
for(g in seq_len(length(PrepList)))
PrepList[[g]]$pars <- pars[[g]]
}
return(PrepList)
}
ReturnPars <- function(PrepList, itemnames, random, lrPars, lr.random = NULL, MG = FALSE){
parnum <- par <- est <- item <- parname <- gnames <- class <-
lbound <- ubound <- prior.type <- prior_1 <- prior_2 <- c()
if(!MG) PrepList <- list(full=PrepList)
for(g in seq_len(length(PrepList))){
tmpgroup <- PrepList[[g]]$pars
for(i in seq_len(length(tmpgroup))){
if(i <= length(itemnames))
item <- c(item, rep(itemnames[i], length(tmpgroup[[i]]@parnum)))
class <- c(class, rep(class(tmpgroup[[i]]), length(tmpgroup[[i]]@parnum)))
parname <- c(parname, names(tmpgroup[[i]]@est))
parnum <- c(parnum, tmpgroup[[i]]@parnum)
par <- c(par, tmpgroup[[i]]@par)
est <- c(est, tmpgroup[[i]]@est)
lbound <- c(lbound, tmpgroup[[i]]@lbound)
ubound <- c(ubound, tmpgroup[[i]]@ubound)
tmp <- sapply(as.character(tmpgroup[[i]]@prior.type),
function(x) switch(x, '1'='norm', '2'='lnorm',
'3'='beta', '4'='expbeta', 'none'))
prior.type <- c(prior.type, tmp)
prior_1 <- c(prior_1, tmpgroup[[i]]@prior_1)
prior_2 <- c(prior_2, tmpgroup[[i]]@prior_2)
}
item <- c(item, rep('GROUP', length(tmpgroup[[i]]@parnum)))
}
for(i in seq_len(length(random))){
parname <- c(parname, names(random[[i]]@est))
parnum <- c(parnum, random[[i]]@parnum)
par <- c(par, random[[i]]@par)
est <- c(est, random[[i]]@est)
lbound <- c(lbound, random[[i]]@lbound)
ubound <- c(ubound, random[[i]]@ubound)
tmp <- sapply(as.character(random[[i]]@prior.type),
function(x) switch(x, '1'='norm', '2'='lnorm',
'3'='beta', '4'='expbeta', 'none'))
prior.type <- c(prior.type, tmp)
prior_1 <- c(prior_1, random[[i]]@prior_1)
prior_2 <- c(prior_2, random[[i]]@prior_2)
class <- c(class, rep('RandomPars', length(random[[i]]@parnum)))
item <- c(item, rep('RANDOM', length(random[[i]]@parnum)))
}
if(length(lrPars)){
parname <- c(parname, names(lrPars@est))
parnum <- c(parnum, lrPars@parnum)
par <- c(par, lrPars@par)
est <- c(est, lrPars@est)
lbound <- c(lbound, lrPars@lbound)
ubound <- c(ubound, lrPars@ubound)
tmp <- sapply(as.character(lrPars@prior.type),
function(x) switch(x, '1'='norm', '2'='lnorm',
'3'='beta', '4'='expbeta', 'none'))
prior.type <- c(prior.type, tmp)
prior_1 <- c(prior_1, lrPars@prior_1)
prior_2 <- c(prior_2, lrPars@prior_2)
class <- c(class, rep('lrPars', length(lrPars@parnum)))
item <- c(item, rep('BETA', length(lrPars@parnum)))
}
for(i in seq_len(length(lr.random))){
parname <- c(parname, names(lr.random[[i]]@est))
parnum <- c(parnum, lr.random[[i]]@parnum)
par <- c(par, lr.random[[i]]@par)
est <- c(est, lr.random[[i]]@est)
lbound <- c(lbound, lr.random[[i]]@lbound)
ubound <- c(ubound, lr.random[[i]]@ubound)
tmp <- sapply(as.character(lr.random[[i]]@prior.type),
function(x) switch(x, '1'='norm', '2'='lnorm',
'3'='beta', '4'='expbeta', 'none'))
prior.type <- c(prior.type, tmp)
prior_1 <- c(prior_1, lr.random[[i]]@prior_1)
prior_2 <- c(prior_2, lr.random[[i]]@prior_2)
class <- c(class, rep('LRRandomPars', length(lr.random[[i]]@parnum)))
item <- c(item, rep('LRRANDOM', length(lr.random[[i]]@parnum)))
}
gnames <- rep(names(PrepList), each = length(est)/length(PrepList))
par[parname %in% c('g', 'u')] <- antilogit(par[parname %in% c('g', 'u')])
lbound[parname %in% c('g', 'u')] <- antilogit(lbound[parname %in% c('g', 'u')])
ubound[parname %in% c('g', 'u')] <- antilogit(ubound[parname %in% c('g', 'u')])
ret <- data.frame(group=gnames, item=item, class=class, name=parname, parnum=parnum, value=par,
lbound=lbound, ubound=ubound, est=est, prior.type=prior.type,
prior_1=prior_1, prior_2=prior_2)
ret
}
UpdatePrepList <- function(PrepList, pars, random, lr.random, lrPars = list(), MG = FALSE){
currentDesign <- ReturnPars(PrepList, PrepList[[1L]]$itemnames, random=random,
lrPars=lrPars, lr.random=lr.random, MG = TRUE)
if(nrow(currentDesign) != nrow(pars))
stop('Rows in supplied and starting value data.frame objects do not match. Were the
data or itemtype input arguments modified?', call.=FALSE)
if(!all(as.matrix(currentDesign[,c('group', 'item', 'class', 'name', 'parnum')]) ==
as.matrix(pars[,c('group', 'item', 'class', 'name', 'parnum')])))
stop('Critical internal parameter labels do not match those returned from pars = \'values\'',
call.=FALSE)
if(!all(sapply(currentDesign, class) == sapply(pars, class)))
stop('pars input does not contain the appropriate classes, which should match pars = \'values\'',
call.=FALSE)
if(!all(unique(pars$prior.type) %in% c('none', 'norm', 'beta', 'lnorm', 'expbeta')))
stop('prior.type input in pars contains invalid prior types', call.=FALSE)
if(!MG) PrepList <- list(PrepList)
pars$value[pars$name %in% c('g', 'u')] <- logit(pars$value[pars$name %in% c('g', 'u')])
pars$lbound[pars$name %in% c('g', 'u')] <- logit(pars$lbound[pars$name %in% c('g', 'u')])
pars$ubound[pars$name %in% c('g', 'u')] <- logit(pars$ubound[pars$name %in% c('g', 'u')])
if(PrepList[[1L]]$nfact > 1L)
PrepList[[1L]]$exploratory <- all(pars$est[pars$name %in% paste0('a', seq_len(PrepList[[1L]]$nfact))])
ind <- 1L
for(g in seq_len(length(PrepList))){
for(i in seq_len(length(PrepList[[g]]$pars))){
for(j in seq_len(length(PrepList[[g]]$pars[[i]]@par))){
PrepList[[g]]$pars[[i]]@par[j] <- pars[ind,'value']
PrepList[[g]]$pars[[i]]@est[j] <- as.logical(pars[ind,'est'])
PrepList[[g]]$pars[[i]]@lbound[j] <- pars[ind,'lbound']
PrepList[[g]]$pars[[i]]@ubound[j] <- pars[ind,'ubound']
tmp <- as.character(pars[ind,'prior.type'])
PrepList[[g]]$pars[[i]]@prior.type[j] <-
switch(tmp, norm=1L, lnorm=2L, beta=3L, expbeta=4L, 0L)
PrepList[[g]]$pars[[i]]@prior_1[j] <- pars[ind,'prior_1']
PrepList[[g]]$pars[[i]]@prior_2[j] <- pars[ind,'prior_2']
ind <- ind + 1L
}
if(is(PrepList[[g]]$pars[[i]], 'graded')){
tmp <- ExtractZetas(PrepList[[g]]$pars[[i]])
if(!all(tmp == sort(tmp, decreasing=TRUE)) || length(unique(tmp)) != length(tmp))
stop('Graded model intercepts for item ', i, ' in group ', g,
' do not descend from highest to lowest. Please fix', call.=FALSE)
}
PrepList[[g]]$pars[[i]]@any.prior <- any(1L:3L %in%
PrepList[[g]]$pars[[i]]@prior.type)
}
}
if(length(random)){
for(i in seq_len(length(random))){
for(j in seq_len(length(random[[i]]@par))){
random[[i]]@par[j] <- pars[ind,'value']
random[[i]]@est[j] <- as.logical(pars[ind,'est'])
random[[i]]@lbound[j] <- pars[ind,'lbound']
random[[i]]@ubound[j] <- pars[ind,'ubound']
ind <- ind + 1L
}
}
attr(PrepList, 'random') <- random
}
if(length(lr.random)){
for(i in seq_len(length(lr.random))){
for(j in seq_len(length(lr.random[[i]]@par))){
lr.random[[i]]@par[j] <- pars[ind,'value']
lr.random[[i]]@est[j] <- as.logical(pars[ind,'est'])
lr.random[[i]]@lbound[j] <- pars[ind,'lbound']
lr.random[[i]]@ubound[j] <- pars[ind,'ubound']
ind <- ind + 1L
}
}
attr(PrepList, 'lr.random') <- lr.random
}
if(!MG) PrepList <- PrepList[[1L]]
return(PrepList)
}
#new gradient and hessian with priors
DerivativePriors <- function(x, grad, hess){
if(any(x@prior.type %in% 1L)){ #norm
ind <- x@prior.type %in% 1L
val <- x@par[ind]
mu <- x@prior_1[ind]
s <- x@prior_2[ind]
g <- -(val - mu)/(s^2)
h <- -1/(s^2)
grad[ind] <- grad[ind] + g
if(length(val) == 1L) hess[ind, ind] <- hess[ind, ind] + h
else diag(hess[ind, ind]) <- diag(hess[ind, ind]) + h
}
if(any(x@prior.type %in% 2L)){ #lnorm
ind <- x@prior.type %in% 2L
val <- x@par[ind]
val <- ifelse(val > 0, val, 1e-10)
lval <- log(val)
mu <- x@prior_1[ind]
s <- x@prior_2[ind]
g <- -(lval - mu)/(val * s^2) - 1/val
h <- 1/(val^2) - 1/(val^2 * s^2) - (lval - mu)/(val^2 * s^2)
grad[ind] <- grad[ind] + g
if(length(val) == 1L) hess[ind, ind] <- hess[ind, ind] + h
else diag(hess[ind, ind]) <- diag(hess[ind, ind]) + h
}
if(any(x@prior.type %in% c(3L, 4L))){ #beta
ind <- x@prior.type %in% c(3L, 4L)
val <- x@par[ind]
ind <- x@prior.type == 4L
val[ind] <- plogis(val[ind])
val <- ifelse(val < 1e-10, 1e-10, val)
val <- ifelse(val > 1-1e-10, 1-1e-10, val)
a <- x@prior_1[ind]
b <- x@prior_2[ind]
g <- (a - 1)/val - (b-1)/(1-val)
h <- -(a - 1)/(val^2) - (b-1) / (1-val)^2
grad[ind] <- grad[ind] + g
if(length(val) == 1L) hess[ind, ind] <- hess[ind, ind] + h
else diag(hess[ind, ind]) <- diag(hess[ind, ind]) + h
}
return(list(grad=grad, hess=hess))
}
#new likelihood with priors
LL.Priors <- function(x, LL){
if(any(x@prior.type %in% 1L)){
ind <- x@prior.type %in% 1L
val <- x@par[ind]
u <- x@prior_1[ind]
s <- x@prior_2[ind]
for(i in seq_len(length(val))){
tmp <- dnorm(val[i], u[i], s[i], log=TRUE)
LL <- LL + ifelse(tmp == -Inf, log(1e-100), tmp)
}
}
if(any(x@prior.type %in% 2L)){
ind <- x@prior.type %in% 2L
val <- x@par[ind]
u <- x@prior_1[ind]
s <- x@prior_2[ind]
for(i in seq_len(length(val))){
if(val[i] > 0)
LL <- LL + dlnorm(val[i], u[i], s[i], log=TRUE)
else LL <- LL + log(1e-100)
}
}
if(any(x@prior.type %in% 3L)){
ind <- x@prior.type %in% 3L
val <- x@par[ind]
a <- x@prior_1[ind]
b <- x@prior_2[ind]
for(i in seq_len(length(val))){
if(val[i] > 0 && val[i] < 1)
LL <- LL + dbeta(val[i], a[i], b[i], log=TRUE)
else LL <- LL + log(1e-100)
}
}
if(any(x@prior.type %in% 4L)){
ind <- x@prior.type %in% 4L
val <- plogis(x@par[ind])
a <- x@prior_1[ind]
b <- x@prior_2[ind]
for(i in seq_len(length(val))){
if(val[i] > 0 && val[i] < 1)
LL <- LL + dbeta(val[i], a[i], b[i], log=TRUE)
else LL <- LL + log(1e-100)
}
}
return(LL)
}
ItemInfo <- function(x, Theta, cosangle, total.info = TRUE){
P <- ProbTrace(x, Theta)
dx <- DerivTheta(x, Theta)
info <- matrix(0, nrow(Theta), ncol(P))
cosanglefull <- matrix(cosangle, nrow(P), length(cosangle), byrow = TRUE)
for(i in seq_len(x@ncat))
dx$grad[[i]] <- matrix(rowSums(dx$grad[[i]] * cosanglefull))
for(i in seq_len(x@ncat))
info[,i] <- (dx$grad[[i]])^2 / P[ ,i]
if(total.info) info <- rowSums(info)
return(info)
}
ItemInfo2 <- function(x, Theta, total.info = TRUE, MD = FALSE, DERIV = NULL, P = NULL){
if(is.null(P)) P <- ProbTrace(x, Theta)
dx <- if(is.null(DERIV)) DerivTheta(x, Theta) else DERIV(x, Theta)
if(MD){
info <- matrix(0, length(Theta), length(Theta))
for(i in seq_len(x@ncat))
info <- info + outer(as.numeric(dx$grad[[i]]), as.numeric(dx$grad[[i]])) / P[ ,i]
} else {
grad <- do.call(cbind, dx$grad)
info <- grad^2 / P
if(total.info) info <- rowSums(info)
}
return(info)
}
nameInfoMatrix <- function(info, correction, L, npars){
#give info meaningful names for wald test
parnames <- names(correction)
tmp <- outer(seq_len(npars), rep(1L, npars))
matind <- matrix(0, ncol(tmp), nrow(tmp))
matind[lower.tri(matind, diag = TRUE)] <- tmp[lower.tri(tmp, diag = TRUE)]
matind <- matind * L
matind[matind == 0 ] <- NA
matind[!is.na(matind)] <- tmp[!is.na(matind)]
shortnames <- c()
for(i in seq_len(length(correction))){
keep <- is.na(matind[ , 1L])
while(all(keep)){
matind <- matind[-1L, -1L, drop = FALSE]
keep <- is.na(matind[ , 1L])
}
tmp <- paste0(parnames[i], paste0('.', matind[!keep, 1L], collapse=''))
shortnames <- c(shortnames, tmp)
matind <- matind[keep, keep, drop = FALSE]
}
colnames(info) <- rownames(info) <- shortnames
return(info)
}
maketabData <- function(stringfulldata, stringtabdata, group, groupNames, nitem, K, itemloc,
Names, itemnames, survey.weights){
tabdata2 <- lapply(strsplit(stringtabdata, split='/'), as.integer)
tabdata2 <- do.call(rbind, tabdata2)
tabdata2[tabdata2 == 99999L] <- NA
tabdata <- matrix(0L, nrow(tabdata2), sum(K))
for(i in seq_len(nitem)){
uniq <- sort(na.omit(unique(tabdata2[,i])))
if(length(uniq) < K[i]) uniq <- 0L:(K[i]-1L)
for(j in seq_len(length(uniq)))
tabdata[,itemloc[i] + j - 1L] <- as.integer(tabdata2[,i] == uniq[j])
}
tabdata[is.na(tabdata)] <- 0L
colnames(tabdata) <- Names
colnames(tabdata2) <- itemnames
groupFreq <- vector('list', length(groupNames))
names(groupFreq) <- groupNames
for(g in seq_len(length(groupNames))){
Freq <- integer(length(stringtabdata))
tmpstringdata <- stringfulldata[group == groupNames[g]]
if(!is.null(survey.weights)){
Freq <- mySapply(seq_len(nrow(tabdata)), function(x, std, tstd, w)
sum(w[stringtabdata[x] == tstd]), std=stringtabdata, tstd=tmpstringdata,
w=survey.weights)
} else {
Freq[stringtabdata %in% tmpstringdata] <- as.integer(table(
match(tmpstringdata, stringtabdata)))
}
groupFreq[[g]] <- Freq
}
ret <- list(tabdata=tabdata, tabdata2=tabdata2, Freq=groupFreq)
ret
}
makeLmats <- function(pars, constrain, random = list(), lrPars = list(), lr.random = list()){
ngroups <- length(pars)
J <- length(pars[[1L]]) - 1L
L <- c()
for(g in seq_len(ngroups))
for(i in seq_len(J+1L))
L <- c(L, pars[[g]][[i]]@est)
for(i in seq_len(length(random)))
L <- c(L, random[[i]]@est)
if(length(lrPars))
L <- c(L, lrPars@est)
for(i in seq_len(length(lr.random)))
L <- c(L, lr.random[[i]]@est)
L <- diag(as.numeric(L))
redun_constr <- rep(FALSE, ncol(L))
for(i in seq_len(length(constrain))){
L[constrain[[i]], constrain[[i]]] <- 1L
for(j in 2L:length(constrain[[i]]))
redun_constr[constrain[[i]][j]] <- TRUE
}
return(list(L=L, redun_constr=redun_constr))
}
updateGrad <- function(g, L) L %*% g
updateHess <- function(h, L) L %*% h %*% L
makeopts <- function(method = 'MHRM', draws = 2000L, calcLL = TRUE, quadpts = NULL,
SE = FALSE, verbose = TRUE, GenRandomPars,
SEtol = .001, grsm.block = NULL, D = 1, TOL = NULL,
rsm.block = NULL, calcNull = FALSE, BFACTOR = FALSE,
technical = list(),
SE.type = 'Oakes', large = NULL, accelerate = 'Ramsay', empiricalhist = FALSE,
optimizer = NULL, solnp_args = list(), nloptr_args = list(), ...)
{
opts <- list()
tnames <- names(technical)
gnames <- c('MAXQUAD', 'NCYCLES', 'BURNIN', 'SEMCYCLES', 'set.seed', 'SEtol', 'symmetric',
'gain', 'warn', 'message', 'customK', 'customPriorFun', 'customTheta', 'MHcand',
'parallel', 'NULL.MODEL', 'theta_lim', 'RANDSTART', 'MHDRAWS', 'removeEmptyRows',
'internal_constraints', 'SEM_window', 'delta', 'MHRM_SE_draws', 'Etable', 'infoAsVcov',
'PLCI', 'plausible.draws', 'storeEtable', 'keep_vcov_PD', 'Norder', 'MCEM_draws')
if(!all(tnames %in% gnames))
stop('The following inputs to technical are invalid: ',
paste0(tnames[!(tnames %in% gnames)], ' '), call.=FALSE)
if((method %in% c('MHRM', 'MIXED', 'SEM')) && SE.type == 'Oakes') SE.type <- 'MHRM'
if(method == 'MCEM' && SE && SE.type != 'complete')
stop('SE.type not currently supported for MCEM method', call.=FALSE)
if(method == 'QMCEM' && SE && SE.type == 'Fisher')
stop('Fisher SE.type not supported for QMCEM method', call.=FALSE)
if((method %in% c('MHRM', 'MIXED', 'SEM')) && !(SE.type %in% c('MHRM', 'FMHRM', 'none')))
stop('SE.type not supported for stochastic method', call.=FALSE)
if(!(method %in% c('MHRM', 'MIXED', 'BL', 'EM', 'QMCEM', 'SEM', 'MCEM')))
stop('method argument not supported', call.=FALSE)
if(!(SE.type %in% c('Richardson', 'forward', 'central', 'crossprod', 'Louis', 'sandwich',
'sandwich.Louis', 'Oakes', 'complete', 'SEM', 'Fisher', 'MHRM', 'FMHRM', 'numerical')))
stop('SE.type argument not supported', call.=FALSE)
if(!(method %in% c('EM', 'QMCEM', 'MCEM'))) accelerate <- 'none'
opts$method = method
if(draws < 1) stop('draws must be greater than 0', call.=FALSE)
opts$draws = draws
opts$theta_lim = technical$theta_lim
opts$calcLL = calcLL
opts$quadpts = quadpts
opts$SE = SE
opts$SE.type = SE.type
opts$verbose = verbose
opts$SEtol = ifelse(is.null(technical$SEtol), .001, technical$SEtol)
opts$grsm.block = grsm.block
opts$rsm.block = rsm.block
opts$calcNull = calcNull
opts$customPriorFun = technical$customPriorFun
opts$dentype <- 'Gaussian'
if(BFACTOR) opts$dentype <- 'bfactor'
if(empiricalhist) opts$dentype <- 'EH'
opts$accelerate = accelerate
opts$Norder <- ifelse(is.null(technical$Norder), 2L, technical$Norder)
opts$delta <- ifelse(is.null(technical$delta), 1e-5, technical$delta)
opts$Etable <- ifelse(is.null(technical$Etable), TRUE, technical$Etable)
opts$plausible.draws <- ifelse(is.null(technical$plausible.draws), 0, technical$plausible.draws)
opts$storeEtable <- ifelse(is.null(technical$storeEtable), FALSE, technical$storeEtable)
if(!is.null(TOL))
if(is.nan(TOL) || is.na(TOL)) opts$calcNull <- FALSE
opts$TOL <- ifelse(is.null(TOL), if(method %in% c('EM', 'QMCEM', 'MCEM')) 1e-4 else
if(method == 'BL') 1e-8 else 1e-3, TOL)
if(SE.type == 'SEM' && SE){
opts$accelerate <- 'none'
if(is.null(TOL)) opts$TOL <- 1e-5
if(is.null(technical$NCYCLES)) technical$NCYCLES <- 1000L
if(method == 'QMCEM')
stop('SEM information matrix not supported with QMCEM estimation', call.=FALSE)
}
if(is.null(technical$symmetric)) technical$symmetric <- TRUE
opts$removeEmptyRows <- if(is.null(technical$removeEmptyRows)) FALSE
else technical$removeEmptyRows
opts$PLCI <- ifelse(is.null(technical$PLCI), FALSE, technical$PLCI)
opts$warn <- if(is.null(technical$warn)) TRUE else technical$warn
opts$message <- if(is.null(technical$message)) TRUE else technical$message
opts$technical <- technical
opts$technical$parallel <- ifelse(is.null(technical$parallel), TRUE, technical$parallel)
opts$MAXQUAD <- ifelse(is.null(technical$MAXQUAD), 20000L, technical$MAXQUAD)
opts$NCYCLES <- ifelse(is.null(technical$NCYCLES), 2000L, technical$NCYCLES)
if(opts$method %in% c('EM', 'QMCEM', 'MCEM'))
opts$NCYCLES <- ifelse(is.null(technical$NCYCLES), 500L, technical$NCYCLES)
opts$BURNIN <- ifelse(is.null(technical$BURNIN), 150L, technical$BURNIN)
opts$SEMCYCLES <- ifelse(is.null(technical$SEMCYCLES), 100L, technical$SEMCYCLES)
opts$SEM_from <- ifelse(is.null(technical$SEM_window), 0, technical$SEM_window[1L])
opts$SEM_to <- ifelse(is.null(technical$SEM_window), 1 - opts$SEtol, technical$SEM_window[2L])
opts$KDRAWS <- ifelse(is.null(technical$KDRAWS), 1L, technical$KDRAWS)
opts$MHDRAWS <- ifelse(is.null(technical$MHDRAWS), 5L, technical$MHDRAWS)
opts$MHRM_SE_draws <- ifelse(is.null(technical$MHRM_SE_draws), 2000L, technical$MHRM_SE_draws)
opts$internal_constraints <- ifelse(is.null(technical$internal_constraints),
TRUE, technical$internal_constraints)
opts$keep_vcov_PD <- ifelse(is.null(technical$keep_vcov_PD), TRUE, technical$keep_vcov_PD)
if(empiricalhist){
if(opts$method != 'EM')
stop('empirical histogram method only applicable when method = \'EM\' ', call.=FALSE)
if(opts$TOL == 1e-4) opts$TOL <- 3e-5
if(is.null(opts$quadpts)) opts$quadpts <- 199L
if(opts$NCYCLES == 500L) opts$NCYCLES <- 2000L
}
if(is.null(opts$theta_lim)) opts$theta_lim <- c(-6,6)
if(method == 'QMCEM' && is.null(opts$quadpts)) opts$quadpts <- 5000L
if((opts$method %in% c('MHRM', 'MIXED') || SE.type == 'MHRM') && !GenRandomPars &&
opts$plausible.draws == 0L)
set.seed(12345L)
if(!is.null(technical$set.seed)) set.seed(technical$set.seed)
opts$gain <- c(0.1, 0.75)
if(!is.null(technical$gain)){
if(length(technical$gain) == 2L && is.numeric(technical$gain))
opts$gain <- technical$gain
}
opts$NULL.MODEL <- ifelse(is.null(technical$NULL.MODEL), FALSE, TRUE)
opts$USEEM <- ifelse(method == 'EM', TRUE, FALSE)
opts$returnPrepList <- FALSE
opts$PrepList <- NULL
if(is.null(optimizer)){
opts$Moptim <- if(method %in% c('EM','BL','QMCEM', 'MCEM')) 'BFGS' else 'NR1'
} else {
opts$Moptim <- optimizer
}
if(method == 'MIXED' && opts$Moptim != 'NR1')
stop('optimizer currently cannot be changed for mixedmirt()', call.=FALSE)
if(opts$Moptim == 'solnp'){
if(is.null(solnp_args$control)) solnp_args$control <- list()
if(is.null(solnp_args$control$trace)) solnp_args$control$trace <- 0
if(!method %in% c('EM', 'QMCEM', 'MCEM'))
stop('solnp only supported for optimization with EM estimation engine',
call.=FALSE)
opts$solnp_args <- solnp_args
} else if(opts$Moptim == 'nloptr'){
if(!method %in% c('EM', 'QMCEM', 'MCEM'))
stop('nloptr only supported for optimization with EM estimation engine',
call.=FALSE)
opts$solnp_args <- nloptr_args
}
if(SE && opts$Moptim %in% c('solnp', 'nloptr')) #TODO
stop('SE computations currently not supported for solnp or nloptr optimizers', call. = FALSE)
if(!is.null(large)){
if(is.logical(large))
if(large) opts$returnPrepList <- TRUE
if(is.list(large)) opts$PrepList <- large
}
if(!is.null(technical$customK)) opts$calcNull <- FALSE
opts$logLik_if_converged <- ifelse(is.null(technical$logLik_if_converged), TRUE,
technical$logLik_if_converged)
opts$info_if_converged <- ifelse(is.null(technical$info_if_converged), TRUE,
technical$info_if_converged)
if(method == 'MCEM'){
opts$accelerate <- 'none'
opts$MCEM_draws <- if(is.null(technical$MCEM_draws))
function(cycles) 500 + (cycles - 1)*2
else technical$MCEM_draws
}
return(opts)
}
reloadPars <- function(longpars, pars, ngroups, J){
return(.Call('reloadPars', longpars, pars, ngroups, J,
attr(pars[[1L]], 'nclasspars')))
}
computeItemtrace <- function(pars, Theta, itemloc, offterm = matrix(0L, 1L, length(itemloc)-1L),
CUSTOM.IND){
itemtrace <- .Call('computeItemTrace', pars, Theta, itemloc, offterm)
if(length(CUSTOM.IND)){
for(i in CUSTOM.IND)
itemtrace[,itemloc[i]:(itemloc[i+1L] - 1L)] <- ProbTrace(pars[[i]], Theta=Theta)
}
return(itemtrace)
}
assignItemtrace <- function(pars, itemtrace, itemloc){
for(i in seq_len(length(pars)-1L))
pars[[i]]@itemtrace <- itemtrace[ ,itemloc[i]:(itemloc[i+1L] - 1L)]
pars[[length(pars)]]@itemtrace <- itemtrace
pars
}
loadESTIMATEinfo <- function(info, ESTIMATE, constrain, warn){
longpars <- ESTIMATE$longpars
pars <- ESTIMATE$pars
ngroups <- length(pars)
J <- length(pars[[1L]]) - 1L
info <- nameInfoMatrix(info=info, correction=ESTIMATE$correction, L=ESTIMATE$L,
npars=length(longpars))
ESTIMATE$info <- info
isna <- is.na(diag(info))
info <- info[!isna, !isna]
acov <- try(solve(info), TRUE)
if(is(acov, 'try-error')){
if(warn)
warning('Could not invert information matrix; model likely is not identified.',
call.=FALSE)
ESTIMATE$fail_invert_info <- TRUE
return(ESTIMATE)
} else ESTIMATE$fail_invert_info <- FALSE
SEtmp <- diag(solve(info))
if(any(SEtmp < 0)){
if(warn)
warning("Negative SEs set to NaN.\n", call.=FALSE)
SEtmp[SEtmp < 0 ] <- NaN
}
SEtmp <- sqrt(SEtmp)
SE <- rep(NA, length(longpars))
SE[ESTIMATE$estindex_unique[!isna]] <- SEtmp
index <- seq_len(length(longpars))
for(i in seq_len(length(constrain)))
SE[index %in% constrain[[i]][-1L]] <- SE[constrain[[i]][1L]]
ind1 <- 1L
for(g in seq_len(ngroups)){
for(i in seq_len(J+1L)){
ind2 <- ind1 + length(pars[[g]][[i]]@par) - 1L
pars[[g]][[i]]@SEpar <- SE[ind1:ind2]
ind1 <- ind2 + 1L
}
}
ESTIMATE$pars <- pars
if(length(ESTIMATE$lrPars))
ESTIMATE$lrPars@SEpar <- SE[ESTIMATE$lrPars@parnum]
return(ESTIMATE)
}
make.randomdesign <- function(random, longdata, covnames, itemdesign, N, LR=FALSE){
ret <- vector('list', length(random))
for(i in seq_len(length(random))){
f <- gsub(" ", "", as.character(random[[i]])[2L])
splt <- strsplit(f, '\\|')[[1L]]
if(any(grepl('\\*', splt[2L]) | grepl('\\+', splt[2L])))
stop('The + and * operators are not supported. Please specify
which effects you want to interact with the : operator, and specify
additional random effects in seperate list elements', call.=FALSE)
gframe <- model.frame(as.formula(paste0('~',splt[2L])), longdata)
sframe <- model.frame(as.formula(paste0('~',splt[1L])), longdata)
levels <- interaction(gframe)
uniq_levels <- unique(levels)
matpar <- diag(1L + ncol(sframe))
estmat <- lower.tri(matpar, diag=TRUE)
ndim <- ncol(matpar)
if(strsplit(f, '+')[[1L]][[1L]] == '-')
estmat[lower.tri(estmat)] <- FALSE
fn <- paste0('COV_', c(splt[2L], colnames(sframe)))
FNCOV <- outer(fn, c(splt[2L], colnames(sframe)), FUN=paste, sep='_')
par <- matpar[lower.tri(matpar, diag=TRUE)]
est <- estmat[lower.tri(estmat, diag=TRUE)]
names(par) <- names(est) <- FNCOV[lower.tri(FNCOV, diag=TRUE)]
drawvals <- matrix(0, length(uniq_levels), ndim,
dimnames=list(uniq_levels, NULL))
mtch <- match(levels, rownames(drawvals))
gdesign <- matrix(1, length(levels), 1L, dimnames = list(NULL, splt[2L]))
if(ncol(sframe) != 0L){
if(grepl('-1+', splt[1L])){
splt[1L] <- strsplit(splt[1L], '-1\\+')[[1]][2]
} else if(grepl('0+', splt[1L]))
splt[1L] <- strsplit(splt[1L], '0\\+')[[1]][2]
gdesign <- cbind(gdesign,
model.matrix(as.formula(paste0('~',splt[1L])), sframe)[,-1L,drop=FALSE])
}
tmp <- matrix(-Inf, ndim, ndim)
diag(tmp) <- 1e-4
lbound <- tmp[lower.tri(tmp, diag=TRUE)]
ret[[i]] <- new('RandomPars',
par=par,
est=est,
SEpar=rep(NaN,length(par)),
ndim=ndim,
lbound=lbound,
ubound=rep(Inf, length(par)),
gframe=gframe,
gdesign=gdesign,
cand.t.var=.5,
any.prior=FALSE,
prior.type=rep(0L, length(par)),
prior_1=rep(NaN,length(par)),
prior_2=rep(NaN,length(par)),
drawvals=drawvals,
mtch=mtch)
}
ret
}
make.lrdesign <- function(df, formula, factorNames, EM=FALSE, TOL){
nfact <- length(factorNames)
if(is.list(formula)){
if(!all(names(formula) %in% factorNames))
stop('List of fixed effect names do not match factor names', call.=FALSE)
estnames <- X <- vector('list', length(formula))
for(i in 1L:length(formula)){
X[[i]] <- model.matrix(formula[[i]], df)
estnames[[i]] <- colnames(X[[i]])
}
X <- do.call(cbind, X)
X <- X[,unique(colnames(X))]
} else {
X <- model.matrix(formula, df)
}
tXX <- t(X) %*% X
inv_tXX <- try(solve(tXX), silent = TRUE)
if(!is.nan(TOL)){
if(is(inv_tXX, 'try-error'))
stop('Latent regression design matrix contains multicollinear terms.', call. = FALSE)
} else inv_tXX <- matrix(0, ncol(tXX), ncol(tXX))
beta <- matrix(0, ncol(X), nfact)
sigma <- matrix(0, nfact, nfact)
diag(sigma) <- 1
if(is.list(formula)){
est <- matrix(FALSE, nrow(beta), ncol(beta))
for(i in 1L:length(formula)){
name <- names(formula)[[i]]
pick <- which(name == factorNames)
est[colnames(X) %in% estnames[[i]], pick] <- TRUE
}
} else est <- matrix(TRUE, nrow(beta), ncol(beta))
est[1L, ] <- FALSE
est <- as.logical(est)
names(est) <- as.character(t(outer(factorNames, colnames(X),
FUN = function(X, Y) paste(X,Y,sep="_"))))
colnames(beta) <- factorNames
rownames(beta) <- colnames(X)
par <- as.numeric(beta)
ret <- new('lrPars',
par=par,
SEpar=rep(NaN,length(par)),
est=est,
beta=beta,
sigma=sigma,
nfact=nfact,
nfixed=ncol(X),
df=df,
X=X,
tXX=tXX,
inv_tXX=inv_tXX,
lbound=rep(-Inf,length(par)),
ubound=rep(Inf,length(par)),
any.prior=FALSE,
prior.type=rep(0L, length(par)),
prior_1=rep(NaN,length(par)),
prior_2=rep(NaN,length(par)),
formula=if(!is.list(formula)) list(formula) else formula,
EM=EM)
ret
}
update.lrPars <- function(df, lrPars){
pick <- df$class == 'lrPars'
df2 <- df[pick, , drop=FALSE]
lrPars@est[] <- df2$est
lrPars@par <- lrPars@beta[] <- df2$value
if(!all(df2$lbound == -Inf))
warning('latent regression parameters cannot be bounded. Ignoring constraint', call.=FALSE)
if(!all(df2$ubound == Inf))
warning('latent regression parameters cannot be bounded. Ignoring constraint', call.=FALSE)
if(!all(df2$prior.type == 'none'))
warning('latent regression parameters do not support prior distribution. Ignoring input.',
call.=FALSE)
lrPars
}
OffTerm <- function(random, J, N){
ret <- numeric(N*J)
for(i in seq_len(length(random))){
tmp <- rowSums(random[[i]]@gdesign*random[[i]]@drawvals[random[[i]]@mtch, ,drop=FALSE])
ret <- ret + tmp
}
return(matrix(ret, N, J))
}
reloadRandom <- function(random, longpars){
for(i in seq_len(length(random))){
parnum <- random[[i]]@parnum
random[[i]]@par <- longpars[min(parnum):max(parnum)]
}
random
}
smooth.cor <- function(x){
eig <- eigen(x)
negvalues <- eig$values <= 0
while (any(negvalues)) {
eig2 <- ifelse(eig$values < 0, 100 * .Machine$double.eps, eig$values)
x <- eig$vectors %*% diag(eig2) %*% t(eig$vectors)
x <- x/sqrt(diag(x) %*% t(diag(x)))
eig <- eigen(x)
negvalues <- eig$values <= .Machine$double.eps
}
x
}
smooth.cov <- function(cov){
ev <- eigen(cov)
v <- ev$values
off <- sum(v) * .01
v[v < 0] <- off
v <- sum(ev$values) * v/sum(v)
v <- ifelse(v < 1e-4, 1e-4, v)
return(ev$vectors %*% diag(v) %*% t(ev$vectors))
}
RMSEA.CI <- function(X2, df, N, ci.lower=.05, ci.upper=.95) {
lower.lambda <- function(lambda) pchisq(X2, df=df, ncp=lambda) - ci.upper
upper.lambda <- function(lambda) pchisq(X2, df=df, ncp=lambda) - ci.lower
lambda.l <- try(uniroot(f=lower.lambda, lower=0, upper=X2)$root, silent=TRUE)
lambda.u <- try(uniroot(f=upper.lambda, lower=0, upper=max(N, X2*5))$root, silent=TRUE)
if(!is(lambda.l, 'try-error')){
RMSEA.lower <- sqrt(lambda.l/(N*df))
} else {
RMSEA.lower <- 0
}
if(!is(lambda.u, 'try-error')){
RMSEA.upper <- sqrt(lambda.u/(N*df))
} else {
RMSEA.upper <- 0
}
return(c(RMSEA.lower, RMSEA.upper))
}
longpars_constrain <- function(longpars, constrain){
for(i in seq_len(length(constrain)))
longpars[constrain[[i]][-1L]] <- longpars[constrain[[i]][1L]]
longpars
}
BL.LL <- function(p, est, longpars, pars, ngroups, J, Theta, PrepList, specific, sitems,
CUSTOM.IND, EHPrior, Data, dentype, itemloc, theta, constrain, lrPars){
longpars[est] <- p
longpars <- longpars_constrain(longpars=longpars, constrain=constrain)
pars2 <- reloadPars(longpars=longpars, pars=pars, ngroups=ngroups, J=J)
gstructgrouppars <- prior <- Prior <- vector('list', ngroups)
full <- length(lrPars) > 0L
if(full){
lrPars@par <- longpars[lrPars@parnum]
lrPars@beta[] <- lrPars@par
lrPars@mus <- lrPars@X %*% lrPars@beta
}
if(dentype == 'EH'){
Prior[[1L]] <- EHPrior[[1L]]
} else if(dentype == 'custom'){
for(g in seq_len(ngroups)){
gp <- pars[[g]][[J+1L]]
Prior[[g]] <- gp@den(gp, Theta)
Prior[[g]] <- Prior[[g]] / sum(Prior[[g]])
}
} else if(dentype == 'discrete'){
for(g in seq_len(ngroups)){
gp <- pars[[g]][[J+1L]]
if(full){
Prior[[g]] <- gp@den(gp, Theta, mus=lrPars@mus)
Prior[[g]] <- Prior[[g]]/rowSums(Prior[[g]])
} else {
Prior[[g]] <- gp@den(gp, Theta)
Prior[[g]] <- Prior[[g]] / sum(Prior[[g]])
}
}
} else {
for(g in seq_len(ngroups)){
gstructgrouppars[[g]] <- ExtractGroupPars(pars2[[g]][[J+1L]])
if(dentype == 'bfactor'){
prior[[g]] <- dnorm(theta, 0, 1)
prior[[g]] <- prior[[g]]/sum(prior[[g]])
Prior[[g]] <- apply(expand.grid(prior[[g]], prior[[g]]), 1L, prod)
next
}
if(full){
Prior[[g]] <- mirt_dmvnorm(Theta[ ,1L:ncol(lrPars@mus),drop=FALSE],
lrPars@mus, gstructgrouppars[[g]]$gcov, quad=TRUE)
Prior[[g]] <- Prior[[g]]/rowSums(Prior[[g]])
} else {
Prior[[g]] <- mirt_dmvnorm(Theta,gstructgrouppars[[g]]$gmeans,
gstructgrouppars[[g]]$gcov)
Prior[[g]] <- Prior[[g]]/sum(Prior[[g]])
}
}
}
LL <- 0
for(g in seq_len(ngroups)){
expected <- Estep.mirt(pars=pars2[[g]],
tabdata=Data$tabdatalong,
freq=if(full) rep(1L, nrow(Prior[[1L]])) else Data$Freq[[g]],
Theta=Theta, prior=Prior[[g]], itemloc=itemloc,
CUSTOM.IND=CUSTOM.IND, full=full, Etable=FALSE)$expected
LL <- LL + sum(Data$Freq[[g]] * log(expected), na.rm = TRUE)
}
LL
}
select_quadpts <- function(nfact) switch(as.character(nfact),
'1'=61, '2'=31, '3'=15, '4'=9, '5'=7, 3)
mirt_rmvnorm <- function(n, mean = rep(0, nrow(sigma)), sigma = diag(length(mean)),
check = FALSE, pre.ev=list())
{
if(!length(pre.ev)){
# Version modified from mvtnorm::rmvnorm, version 0.9-9996, 19-April, 2014.
if(check){
if (!isSymmetric(sigma, tol = sqrt(.Machine$double.eps), check.attributes = FALSE))
stop("sigma must be a symmetric matrix", call.=FALSE)
if (length(mean) != nrow(sigma))
stop("mean and sigma have non-conforming size", call.=FALSE)
}
ev <- eigen(sigma, symmetric = TRUE)
NCOL <- ncol(sigma)
if(check)
if (!all(ev$values >= -sqrt(.Machine$double.eps) * abs(ev$values[1])))
warning("sigma is numerically not positive definite", call.=FALSE)
} else {
ev <- pre.ev
NCOL <- length(ev$values)
}
retval <- ev$vectors %*% diag(sqrt(ev$values), NCOL) %*% t(ev$vectors)
retval <- matrix(rnorm(n * NCOL), nrow = n) %*% retval
retval <- sweep(retval, 2L, mean, "+")
colnames(retval) <- names(mean)
retval
}
mirt_dmvnorm <- function(x, mean, sigma, log = FALSE, quad = FALSE, stable = TRUE, ...)
{
if(quad && is.matrix(mean)){
isigma <- solve(sigma)
distval <- matrix(0, nrow(mean), nrow(x))
for(i in seq_len(nrow(mean))){
centered <- t(t(x) - mean[i,])
distval[i, ] <- rowSums((centered %*% isigma) * centered)
}
} else {
if(is.matrix(mean)){
centered <- x - mean
distval <- rowSums((centered %*% solve(sigma)) * centered)
} else {
distval <- mahalanobis(x, center = mean, cov = sigma)
}
}
logdet <- sum(log(eigen(sigma, symmetric=TRUE,
only.values=TRUE)$values))
logretval <- -(ncol(x)*log(2*pi) + logdet + distval)/2
if(stable)
logretval <- ifelse(logretval < -690.7755, -690.7755, logretval)
if(log) return(logretval)
exp(logretval)
}
# prior for latent class analysis
lca_prior <- function(Theta, Etable){
TP <- nrow(Theta)
if ( is.null(Etable) ){
prior <- rep( 1/TP , TP )
} else {
prior <- rowSums(Etable)
}
prior <- prior / sum(prior)
return(prior)
}
makeObstables <- function(dat, K, which.items){
ret <- vector('list', ncol(dat))
sumscore <- rowSums(dat)
for(i in seq_len(length(ret))){
if(!(i %in% which.items)) next
ret[[i]] <- matrix(0, sum(K-1L)+1L, K[i])
colnames(ret[[i]]) <- paste0(1L:K[i]-1L)
rownames(ret[[i]]) <- paste0(1L:nrow(ret[[i]])-1L)
split <- by(sumscore, dat[,i], table)
for(j in seq_len(length(split))){
m <- match(names(split[[j]]), rownames(ret[[i]]))
ret[[i]][m,j] <- split[[j]]
}
ret[[i]] <- ret[[i]][-c(1L, nrow(ret[[i]])), ]
}
ret
}
collapseCells <- function(O, E, mincell = 1){
for(i in seq_len(length(O))){
On <- O[[i]]
En <- E[[i]]
if(is.null(En)) next
drop <- which(rowSums(is.na(En)) > 0)
En[is.na(En)] <- 0
#collapse known upper and lower sparce cells
if(length(drop)){
up <- drop[1L]:drop[length(drop)/2]
low <- drop[length(drop)/2 + 1L]:drop[length(drop)]
En[max(up)+1, ] <- colSums(En[c(up, max(up)+1), , drop = FALSE])
On[max(up)+1, ] <- colSums(On[c(up, max(up)+1), , drop = FALSE])
En[min(low)-1, ] <- colSums(En[c(low, min(low)-1), , drop = FALSE])
On[min(low)-1, ] <- colSums(On[c(low, min(low)-1), , drop = FALSE])
En[c(up, low), ] <- On[c(up, low), ] <- NA
En <- na.omit(En)
On <- na.omit(On)
}
#drop 0's and 1's
drop <- rowSums(On) == 0L
On <- On[!drop,]
En <- En[!drop,]
L <- En < mincell
drop <- c()
for(j in seq_len(nrow(On)-1L)){
ss <- sum(On[j,])
if(ss == 1L){
drop <- c(drop, j)
On[j+1L, ] <- On[j+1L, ] + On[j, ]
En[j+1L, ] <- En[j+1L, ] + En[j, ]
}
}
if(length(drop)){
On <- On[-drop,]
En <- En[-drop,]
}
ss <- sum(On[nrow(On),])
if(ss == 1L){
On[nrow(On)-1L, ] <- On[nrow(On)-1L, ] + On[nrow(On), ]
En[nrow(On)-1L, ] <- En[nrow(On)-1L, ] + En[nrow(On), ]
On <- On[-nrow(On),]; En <- En[-nrow(En),]
}
#collapse accross as much as possible
if(ncol(En) > 2L){
for(j in seq_len(nrow(En))){
if(!any(L[j,])) next
tmp <- En[j, ]
tmp2 <- On[j, ]
while(length(tmp) > 2L){
m <- min(tmp)
whc <- max(which(m == tmp))
if(whc == 1L){
tmp[2L] <- tmp[2L] + tmp[1L]
tmp2[2L] <- tmp2[2L] + tmp2[1L]
} else if(whc == length(tmp)){
tmp[length(tmp)-1L] <- tmp[length(tmp)-1L] + tmp[length(tmp)]
tmp2[length(tmp2)-1L] <- tmp2[length(tmp2)-1L] + tmp2[length(tmp2)]
} else {
left <- min(tmp[whc-1L], tmp[whc+1L]) == c(tmp[whc-1L], tmp[whc+1L])[1L]
pick <- if(left) whc-1L else whc+1L
tmp[pick] <- tmp[pick] + tmp[whc]
tmp2[pick] <- tmp2[pick] + tmp2[whc]
}
tmp[whc] <- tmp2[whc] <- NA
tmp <- na.omit(tmp); tmp2 <- na.omit(tmp2)
if(all(tmp >= mincell)) break
}
tmp <- c(tmp, rep(NA, ncol(En)-length(tmp)))
tmp2 <- c(tmp2, numeric(ncol(En)-length(tmp2)))
En[j, ] <- tmp
On[j, ] <- tmp2
}
}
En[is.na(En)] <- 0
# collapse right columns if they are too rare
if(ncol(En) > 2L){
while(TRUE){
pick <- colSums(En) < mincell * ceiling(nrow(En) * .1)
if(!pick[length(pick)] || ncol(En) == 2L) break
if(pick[length(pick)]){
On[ ,length(pick)-1L] <- On[ ,length(pick)-1L] + On[ ,length(pick)]
En[ ,length(pick)-1L] <- En[ ,length(pick)-1L] + En[ ,length(pick)]
On <- On[ ,-length(pick)]; En <- En[ ,-length(pick)]
}
}
}
dropcol <- logical(ncol(En))
#drop all other columns if they are very rare
for(j in ncol(En):2L){
tmp <- sum(En[,j] > 0) / nrow(En)
if(tmp < .05){
dropcol[j] <- TRUE
En[,j-1L] <- En[,j-1L] + En[,j]
On[,j-1L] <- On[,j-1L] + On[,j]
}
}
En <- En[,!dropcol]; On <- On[,!dropcol]
#merge across
L <- En < mincell & En != 0
while(any(L, na.rm = TRUE)){
if(!is.matrix(L)) break
whc <- min(which(rowSums(L) > 0L))
if(whc == 1L){
En[2L,] <- En[2L, ] + En[1L,]
On[2L,] <- On[2L, ] + On[1L,]
En <- En[-1L,]; On <- On[-1L,]
} else if(whc == nrow(En)){
En[nrow(En)-1L,] <- En[nrow(En)-1L, ] + En[nrow(En),]
On[nrow(On)-1L,] <- On[nrow(On)-1L, ] + On[nrow(On),]
En <- En[-nrow(En),]; On <- On[-nrow(On),]
} else {
ss <- c(sum(On[whc-1L,]), sum(On[whc+1L,]))
up <- (min(ss) == ss)[1L]
pick <- if(up) whc-1L else whc+1L
En[pick,] <- En[pick, ] + En[whc,]
On[pick,] <- On[pick, ] + On[whc,]
En <- En[-whc,]; On <- On[-whc,]
}
L <- En < mincell & En != 0
}
En[En == 0] <- NA
E[[i]] <- En
O[[i]] <- On
}
return(list(O=O, E=E))
}
MGC2SC <- function(x, which){
tmp <- x@ParObjects$pars[[which]]
tmp@Model$lrPars <- x@ParObjects$lrPars
ind <- 1L
for(i in seq_len(x@Data$nitems)){
tmp@ParObjects$pars[[i]]@parnum[] <- seq(ind, ind + length(tmp@ParObjects$pars[[i]]@parnum) - 1L)
ind <- ind + length(tmp@ParObjects$pars[[i]]@parnum)
}
tmp@Data <- x@Data
tmp@Data$data <- tmp@Data$data[tmp@Data$group == tmp@Data$groupName[which], , drop=FALSE]
tmp@Data$Freq[[1L]] <- tmp@Data$Freq[[which]]
tmp@Data$fulldata[[1L]] <- x@Data$fulldata[[which]]
tmp
}
#' Compute numerical derivatives
#'
#' Compute numerical derivatives using forward/backword difference,
#' central difference, or Richardson extropolation.
#'
#' @param par a vector of parameters
#' @param f the objective function being evaluated
#' @param ... additional arguments to be passed to \code{f} and the \code{numDeriv} package when the
#' Richardson type is used
#' @param delta the term used to perturb the \code{f} function. Default is 1e-5
#' @param gradient logical; compute the gradient terms? If FALSE then the Hessian is computed instead
#' @param type type of difference to compute. Can be either \code{'forward'} for the forward difference,
#' \code{'central'} for the central difference, or \code{'Richardson'} for the Richardson extropolation.
#' Backword difference is acheived by supplying a negative \code{delta} value
#' @export numerical_deriv
#' @author Phil Chalmers \email{rphilip.chalmers@@gmail.com}
#' @keywords numerical derivatives
#'
#' @examples
#'
#' \dontrun{
#' f <- function(x) 3*x[1]^3 - 4*x[2]^2
#' par <- c(3,8)
#'
#' # grad = 9 * x^2 , -8 * y
#' (actual <- c(9 * par[1]^2, -8 * par[2]))
#' numerical_deriv(par, f, type = 'forward')
#' numerical_deriv(par, f, type = 'central')
#' numerical_deriv(par, f, type = 'Richardson')
#'
#' # hessian = h11 -> 18 * x, h22 -> -8, h12 -> h21 -> 0
#' (actual <- matrix(c(18 * par[1], 0, 0, -8), 2, 2))
#' numerical_deriv(par, f, type = 'forward', gradient = FALSE)
#' numerical_deriv(par, f, type = 'central', gradient = FALSE)
#' numerical_deriv(par, f, type = 'Richardson', gradient = FALSE)
#'
#' }
numerical_deriv <- function(par, f, ..., delta = 1e-5, gradient = TRUE, type = 'forward'){
forward_difference <- function(par, f, delta, ...){
dots <- list(...)
np <- length(par)
g <- numeric(np)
if(is.null(dots$ObJeCtIvE)) fx <- f(par, ...) else fx <- dots$ObJeCtIvE
for(i in seq_len(np)){
p <- par
p[i] <- p[i] + delta
g[i] <- (f(p, ...) - fx) / delta
}
g
}
forward_difference2 <- function(par, f, delta, ...){
dots <- list(...)
np <- length(par)
hess <- matrix(0, np, np)
if(is.null(dots$ObJeCtIvE)) fx <- f(par, ...) else fx <- dots$ObJeCtIvE
fx1 <- numeric(np)
for(i in seq_len(np)){
tmp <- par
tmp[i] <- tmp[i] + delta
fx1[i] <- f(tmp, ...)
}
for(i in seq_len(np)){
for(j in i:np){
fx1x2 <- par
fx1x2[i] <- fx1x2[i] + delta
fx1x2[j] <- fx1x2[j] + delta
hess[i,j] <- hess[j, i] <- (f(fx1x2, ...) - fx1[i] - fx1[j] + fx) / (delta^2)
}
}
hess
}
central_difference <- function(par, f, delta, ...){
np <- length(par)
g <- numeric(np)
for(i in seq_len(np)){
p1 <- p2 <- par
p1[i] <- p1[i] + delta
p2[i] <- p2[i] - delta
g[i] <- (f(p1, ...) - f(p2, ...)) / (2 * delta)
}
g
}
central_difference2 <- function(par, f, delta, ...){
np <- length(par)
hess <- matrix(0, np, np)
fx <- f(par, ...)
for(i in seq_len(np)){
for(j in i:np){
if(i == j){
p1 <- p2 <- par
p1[i] <- p1[i] + delta; s2 <- f(p1, ...)
p1[i] <- p1[i] + delta; s1 <- f(p1, ...)
p2[i] <- p2[i] - delta; s3 <- f(p2, ...)
p2[i] <- p2[i] - delta; s4 <- f(p2, ...)
hess[i, i] <- (-s1 + 16 * s2 - 30 * fx + 16 * s3 - s4) / (12 * delta^2)
} else {
p <- par
p[i] <- p[i] + delta; p[j] <- p[j] + delta; s1 <- f(p, ...)
p[j] <- p[j] - 2*delta; s2 <- f(p, ...)
p[i] <- p[i] - 2*delta; s4 <- f(p, ...)
p[j] <- p[j] + 2*delta; s3 <- f(p, ...)
hess[i,j] <- hess[j,i] <- (s1 - s2 - s3 + s4) / (4 * delta^2)
}
}
}
hess
}
if(!length(par)){
if(gradient) return(numeric())
else return(matrix(numeric()))
}
if(type == 'central'){
ret <- if(gradient) central_difference(par=par, f=f, delta=delta/2, ...)
else central_difference2(par=par, f=f, delta=delta, ...)
} else if(type == 'forward'){
ret <- if(gradient) forward_difference(par=par, f=f, delta=delta, ...)
else forward_difference2(par=par, f=f, delta=delta, ...)
} else if(type == 'Richardson'){
ret <- if(gradient) numDeriv::grad(f, par, ...)
else numDeriv::hessian(f, par, ...)
}
ret
}
computeNullModel <- function(data, itemtype, key, group=NULL){
if(is.null(itemtype)) itemtype <- rep('graded', ncol(data))
itemtype[itemtype == 'Rasch'] <- 'gpcm'
if(!is.null(group)){
null.mod <- multipleGroup(data, 1L, itemtype=itemtype, group=group, verbose=FALSE,
key=key, technical=list(NULL.MODEL=TRUE))
} else {
null.mod <- mirt(data, 1L, itemtype=itemtype, verbose=FALSE, key=key,
technical=list(NULL.MODEL=TRUE))
}
null.mod
}
loadSplineParsItem <- function(x, Theta){
sargs <- x@sargs
Theta_prime <- if(x@stype == 'bs'){
splines::bs(Theta, df=sargs$df, knots=sargs$knots,
degree=sargs$degree, intercept=sargs$intercept)
} else if(x@stype == 'ns'){
splines::ns(Theta, df=sargs$df, knots=sargs$knots,
intercept=sargs$intercept)
}
class(Theta_prime) <- 'matrix'
x@Theta_prime <- Theta_prime
x
}
loadSplinePars <- function(pars, Theta, MG = TRUE){
fn <- function(pars, Theta){
cls <- sapply(pars, class)
pick <- which(cls == 'spline')
if(length(pick)){
for(i in pick)
pars[[i]] <- loadSplineParsItem(pars[[i]], Theta)
}
return(pars)
}
if(MG){
for(g in seq_len(length(pars))){
pars[[g]] <- fn(pars[[g]], Theta)
}
} else {
pars <- fn(pars, Theta)
}
return(pars)
}
latentRegression_obj <- function(data, covdata, formula, empiricalhist, method){
if(!is.null(covdata) && !is.null(formula)){
if(empiricalhist)
stop('Empirical histogram method not supported with covariates', call.=FALSE)
if(!is.data.frame(covdata))
stop('covdata must be a data.frame object', call.=FALSE)
if(nrow(covdata) != nrow(data))
stop('number of rows in covdata do not match number of rows in data', call.=FALSE)
if(!(method %in% c('EM', 'QMCEM')))
stop('method must be from the EM estimation family', call.=FALSE)
tmp <- apply(covdata, 1, function(x) sum(is.na(x)) > 0)
if(any(tmp)){
message('removing rows with NAs in covdata')
covdata <- covdata[-tmp, ]
data <- data[-tmp, ]
}
latent.regression <- list(df=covdata, formula=formula, EM=TRUE)
} else latent.regression <- NULL
latent.regression
}
# borrowed and modified from emdbook package, March 1 2017
mixX2 <- function (p, df = 1, mix = 0.5, lower.tail = TRUE)
{
df <- rep(df, length.out = length(p))
mix <- rep(mix, length.out = length(p))
c1 <- ifelse(df == 1, if (lower.tail) 1
else 0, pchisq(p, df - 1, lower.tail = lower.tail))
c2 <- pchisq(p, df, lower.tail = lower.tail)
r <- mix * c1 + (1 - mix) * c2
r
}
get_deriv_coefs <- function(order, deriv = 1L){
if(deriv == 1L){
ret <- switch(as.character(order),
"1" = c(-1, 1),
"2" = c(-1/2, 1/2))
}
ret
}
cfi <- function(X2, X2.null, df, df.null){
ret <- 1 - (X2 - df) / (X2.null - df.null)
if(ret > 1) ret <- 1
if(ret < 0) ret <- 0
ret
}
tli <- function(X2, X2.null, df, df.null)
(X2.null/df.null - X2/df) / (X2.null/df.null - 1)
rmsea <- function(X2, df, N){
ret <- ifelse((X2 - df) > 0,
sqrt(X2 - df) / sqrt(df * (N-1)), 0)
ret <- ifelse(is.na(ret), NaN, ret)
ret
}
controlCandVar <- function(PA, cand, min = .1, max = .6){
if(PA > max) cand <- cand * 1.05
else if(PA < min) cand <- cand * 0.9
if(cand < .001) cand <- .001
cand
}
QMC_quad <- function(npts, nfact, lim, leap=409, norm=FALSE){
qnorm(sfsmisc::QUnif(npts, min=0, max=1, p=nfact, leap=leap))
}
MC_quad <- function(npts, nfact, lim)
qnorm(matrix(runif(n=npts * nfact, min = lim[1L], max = lim[2]), npts, nfact))
respSample <- function(P) .Call("respSample", P)
missingMsg <- function(string)
stop(paste0('\'', string, '\' argument is missing.'), call.=FALSE)
.mirtClusterEnv <- new.env(parent=emptyenv())
.mirtClusterEnv$ncores <- 1L
myApply <- function(X, MARGIN, FUN, ...){
if(.mirtClusterEnv$ncores > 1L){
return(t(parallel::parApply(cl=.mirtClusterEnv$MIRTCLUSTER, X=X, MARGIN=MARGIN, FUN=FUN, ...)))
} else {
return(t(apply(X=X, MARGIN=MARGIN, FUN=FUN, ...)))
}
}
myLapply <- function(X, FUN, ...){
if(.mirtClusterEnv$ncores > 1L){
return(parallel::parLapply(cl=.mirtClusterEnv$MIRTCLUSTER, X=X, fun=FUN, ...))
} else {
return(lapply(X=X, FUN=FUN, ...))
}
}
mySapply <- function(X, FUN, ...){
if(.mirtClusterEnv$ncores > 1L){
return(t(parallel::parSapply(cl=.mirtClusterEnv$MIRTCLUSTER, X=X, FUN=FUN, ...)))
} else {
return(t(sapply(X=X, FUN=FUN, ...)))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.