Skip to content

Commit df644e0

Browse files
philchalmersphilchalmers
committed
support Oakes method for everything
1 parent 8fd62f8 commit df644e0

File tree

7 files changed

+88
-28
lines changed

7 files changed

+88
-28
lines changed

‎NEWS.md

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,8 @@
11
# Changes in mirt 1.43
22

3+
- Added support for latent regression ACOV/SE estimation with Oakes method
4+
in `mirt()`
5+
36
- Related to both points below, general MLTM (Embretson, 1984) added when
47
itemtype is specified as `PC1PL` and an `itemdesign` set is used, where
58
formula must include the name of the factor in the formula expressions. See

‎R/03-estimation.R

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -625,9 +625,9 @@ ESTIMATION <- function(data, model, group, itemtype = NULL, guess = 0, upper = 1
625625
if(opts$method %in% c('EM', 'BL', 'QMCEM', 'MCEM')){
626626
logLik <- G2 <- SElogLik <- 0
627627
if(length(lrPars)){
628-
if(opts$SE && !(opts$SE.type %in% c('complete', 'forward', 'central', 'Richardson')))
629-
stop('Information matrix method for latent regression estimates not supported',
630-
call.=FALSE)
628+
# if(opts$SE && !(opts$SE.type %in% c('complete', 'forward', 'central', 'Richardson')))
629+
# stop('Information matrix method for latent regression estimates not supported',
630+
# call.=FALSE)
631631
opts$full <- TRUE
632632
} else opts$full <- FALSE
633633
temp <- matrix(0L,nrow=nitems,ncol=nspec)

‎R/08-SE.methods.R

Lines changed: 23 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -223,7 +223,7 @@ SE.Oakes <- function(pick, pars, L, constrain, est, shortpars, longpars,
223223
rlist, full, Data, specific, itemloc, CUSTOM.IND,
224224
delta, prior, Prior, Priorbetween, nfact, mixtype,
225225
PrepList, ANY.PRIOR, DERIV, SLOW.IND, Norder, omp_threads,
226-
zero_g = NULL){
226+
lrPars, zero_g = NULL){
227227
r <- 1L
228228
Richardson <- if(Norder > 2L) TRUE else FALSE
229229
if(Richardson){
@@ -249,7 +249,12 @@ SE.Oakes <- function(pick, pars, L, constrain, est, shortpars, longpars,
249249
longpars <- longpars_constrain(longpars, constrain)
250250
pars <- reloadPars(longpars=longpars, pars=pars,
251251
ngroups=ngroups, J=J)
252-
tmp <- updatePrior(pars=pars, gTheta=gTheta,
252+
if(length(lrPars)){
253+
lrPars@par <- longpars[lrPars@parnum]
254+
lrPars@beta[] <- matrix(lrPars@par, lrPars@nfixed, lrPars@nfact)
255+
lrPars@mus <- lrPars@X %*% lrPars@beta
256+
}
257+
tmp <- updatePrior(pars=pars, gTheta=gTheta, lrPars=lrPars,
253258
list=list, ngroups=ngroups, nfact=nfact,
254259
J=J, dentype=dentype, sitems=sitems, cycles=100L,
255260
rlist=rlist, full=full, MC=list$method == 'QMCEM')
@@ -263,6 +268,11 @@ SE.Oakes <- function(pick, pars, L, constrain, est, shortpars, longpars,
263268
longpars <- longpars_old
264269
pars <- reloadPars(longpars=longpars, pars=pars,
265270
ngroups=ngroups, J=J)
271+
if(length(lrPars)){
272+
lrPars@par <- longpars[lrPars@parnum]
273+
lrPars@beta[] <- matrix(lrPars@par, lrPars@nfixed, lrPars@nfact)
274+
lrPars@mus <- lrPars@X %*% lrPars@beta
275+
}
266276
if(pars[[1L]][[J + 1L]]@itemclass == -1L){
267277
for(g in seq_len(length(pars))){
268278
gp <- pars[[g]][[J + 1L]]
@@ -304,6 +314,17 @@ SE.Oakes <- function(pick, pars, L, constrain, est, shortpars, longpars,
304314
deriv <- Deriv.mix(mixtype)
305315
g[mixtype$parnum] <- deriv$grad
306316
}
317+
if(length(lrPars)){
318+
for(group in seq_len(ngroups)){
319+
gp <- ExtractGroupPars(pars[[group]][[J+1L]])
320+
tmp <- Mstep.LR(Theta=gTheta[[group]], CUSTOM.IND=CUSTOM.IND, pars=pars[[group]],
321+
itemloc=itemloc, fulldata=Data$fulldata[[1L]], prior=Prior[[group]],
322+
lrPars=lrPars, retscores=TRUE)
323+
deriv <- Deriv(lrPars, cov=gp$gcov, theta=tmp)
324+
deriv$grad * lrPars@est
325+
g[lrPars@parnum] <- as.vector(deriv$grad * lrPars@est)
326+
}
327+
}
307328
tmp <- g %*% L
308329
if(pick == 0L) return(tmp[est])
309330
grad[row, ] <- tmp[est]

‎R/EMstep.group.R

Lines changed: 14 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -102,8 +102,11 @@ EM.group <- function(pars, constrain, Ls, Data, PrepList, list, Theta, DERIV, so
102102
for(g in seq_len(ngroups)){
103103
for(j in seq_len(J+1L))
104104
est <- c(est, pars[[g]][[j]]@est)
105-
if(length(lrPars))
106-
est <- c(est, rep(FALSE, length(lrPars@est)))
105+
if(length(lrPars)){
106+
tmp <- rep(FALSE, length(lrPars@est))
107+
names(tmp) <- names(lrPars@est)
108+
est <- c(est, tmp)
109+
}
107110
}
108111
for(i in seq_len(length(constrain)))
109112
est[constrain[[i]][-1L]] <- FALSE
@@ -481,6 +484,7 @@ EM.group <- function(pars, constrain, Ls, Data, PrepList, list, Theta, DERIV, so
481484
for(i in 0L:(ncol(deriv$grad)-1L))
482485
h[lrPars@parnum[1L:nrow(deriv$grad) + nrow(deriv$grad)*i],
483486
lrPars@parnum[1L:nrow(deriv$grad) + nrow(deriv$grad)*i]] <- deriv$hess
487+
est[lrPars@parnum] <- lrPars@est
484488
}
485489
}
486490
if(dentype == 'mixture'){
@@ -498,9 +502,10 @@ EM.group <- function(pars, constrain, Ls, Data, PrepList, list, Theta, DERIV, so
498502
} else mixtype <- NULL
499503
hess <- updateHess(h=h, L=L)
500504
hess <- as.matrix(hess[estpars & !redun_constr, estpars & !redun_constr])
501-
if(list$SE.type %in% c('Oakes', 'sandwich') && length(lrPars) && list$SE){
502-
warning('Oakes method not supported for models with latent regression effects', call.=FALSE)
503-
} else if(list$SE.type %in% c('Oakes', 'sandwich') && list$SE){
505+
# if(list$SE.type %in% c('Oakes', 'sandwich') && length(lrPars) && list$SE){
506+
# warning('Oakes method not supported for models with latent regression effects', call.=FALSE)
507+
# } else
508+
if(list$SE.type %in% c('Oakes', 'sandwich') && list$SE){
504509
complete_info <- hess
505510
shortpars <- longpars[estpars & !redun_constr]
506511
tmp <- updatePrior(pars=pars, gTheta=gTheta,
@@ -513,7 +518,7 @@ EM.group <- function(pars, constrain, Ls, Data, PrepList, list, Theta, DERIV, so
513518
pars=pars, L=L, constrain=constrain, delta=list$delta,
514519
est=est, shortpars=shortpars, longpars=longpars,
515520
gTheta=gTheta, list=list, ngroups=ngroups, J=J,
516-
dentype=dentype, sitems=sitems, nfact=nfact,
521+
dentype=dentype, sitems=sitems, nfact=nfact, lrPars=lrPars,
517522
rlist=rlist, full=full, Data=Data, mixtype=mixtype,
518523
specific=specific, itemloc=itemloc, CUSTOM.IND=CUSTOM.IND,
519524
prior=prior, Priorbetween=Priorbetween, Prior=Prior,
@@ -523,7 +528,7 @@ EM.group <- function(pars, constrain, Ls, Data, PrepList, list, Theta, DERIV, so
523528
zero_g <- SE.Oakes(0L, pars=pars, L=L, constrain=constrain, delta=0,
524529
est=est, shortpars=shortpars, longpars=longpars,
525530
gTheta=gTheta, list=list, ngroups=ngroups, J=J,
526-
dentype=dentype, sitems=sitems, nfact=nfact,
531+
dentype=dentype, sitems=sitems, nfact=nfact, lrPars=lrPars,
527532
rlist=rlist, full=full, Data=Data, mixtype=mixtype,
528533
specific=specific, itemloc=itemloc, CUSTOM.IND=CUSTOM.IND,
529534
prior=prior, Priorbetween=Priorbetween, Prior=Prior,
@@ -544,6 +549,8 @@ EM.group <- function(pars, constrain, Ls, Data, PrepList, list, Theta, DERIV, so
544549
pars <- reloadPars(longpars=longpars, pars=pars,
545550
ngroups=ngroups, J=J)
546551
is.latent <- grepl('MEAN_', names(shortpars)) | grepl('COV_', names(shortpars))
552+
if(length(lrPars))
553+
is.latent[names(shortpars) %in% names(lrPars@est)] <- TRUE
547554
missing_info[is.latent, is.latent] <- 0
548555
hess <- complete_info + missing_info
549556
}

‎R/mirt.R

Lines changed: 22 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1191,23 +1191,37 @@
11911191
#'
11921192
#' # covariates
11931193
#' X1 <- rnorm(N); X2 <- rnorm(N)
1194-
#' covdata <- data.frame(X1, X2)
1194+
#' covdata <- data.frame(X1, X2, X3 = rnorm(N))
11951195
#' Theta <- matrix(0.5 * X1 + -1 * X2 + rnorm(N, sd = 0.5))
11961196
#'
11971197
#' # items and response data
11981198
#' a <- matrix(1, 20); d <- matrix(rnorm(20))
11991199
#' dat <- simdata(a, d, 1000, itemtype = '2PL', Theta=Theta)
12001200
#'
12011201
#' # unconditional Rasch model
1202-
#' mod0 <- mirt(dat, 1, 'Rasch')
1202+
#' mod0 <- mirt(dat, 1, 'Rasch', SE=TRUE)
1203+
#' coef(mod0, printSE=TRUE)
12031204
#'
1204-
#' # conditional model using X1 and X2 as predictors of Theta
1205-
#' mod1 <- mirt(dat, 1, 'Rasch', covdata=covdata, formula = ~ X1 + X2)
1205+
#' # conditional model using X1, X2, and X3 (bad) as predictors of Theta
1206+
#' mod1 <- mirt(dat, 1, 'Rasch', covdata=covdata, formula = ~ X1 + X2 + X3, SE=TRUE)
1207+
#' coef(mod1, printSE=TRUE)
12061208
#' coef(mod1, simplify=TRUE)
1207-
#' anova(mod0, mod1)
1208-
#'
1209-
#' # bootstrapped confidence intervals
1210-
#' boot.mirt(mod1, R=5)
1209+
#' anova(mod0, mod1) # jointly significant predictors of theta
1210+
#'
1211+
#' # large sample z-ratios and p-values (if one cares)
1212+
#' cfs <- coef(mod1, printSE=TRUE)
1213+
#' (z <- cfs$lr.betas[[1]] / cfs$lr.betas[[2]])
1214+
#' round(pnorm(abs(z[,1]), lower.tail=FALSE)*2, 3)
1215+
#'
1216+
#' # drop predictor for nested comparison
1217+
#' mod1b <- mirt(dat, 1, 'Rasch', covdata=covdata, formula = ~ X1 + X2)
1218+
#' anova(mod1b, mod1)
1219+
#'
1220+
#' # compare to mixedmirt() version of the same model
1221+
#' mod1.mixed <- mixedmirt(dat, 1, itemtype='Rasch',
1222+
#' covdata=covdata, lr.fixed = ~ X1 + X2 + X3, SE=TRUE)
1223+
#' coef(mod1.mixed)
1224+
#' coef(mod1.mixed, printSE=TRUE)
12111225
#'
12121226
#' # draw plausible values for secondary analyses
12131227
#' pv <- fscores(mod1, plausible.draws = 10)

‎R/utils.R

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2000,6 +2000,7 @@ make.lrdesign <- function(df, formula, factorNames, EM=FALSE, TOL){
20002000
colnames(beta) <- factorNames
20012001
rownames(beta) <- colnames(X)
20022002
par <- as.numeric(beta)
2003+
names(par) <- names(est)
20032004
ret <- new('lrPars',
20042005
par=par,
20052006
SEpar=rep(NaN,length(par)),

‎man/mirt.Rd

Lines changed: 22 additions & 8 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

0 commit comments

Comments
 (0)