Mixed effects modeling for MIRT models

mixedmirt fits MIRT models using FIML estimation to dichotomous and polytomous IRT models conditional on fixed and random effect of person and item level covariates. This can also be understood as 'explanatory IRT' if only fixed effects are modeled, or multilevel/mixed IRT if random and fixed effects are included. The method uses the MH-RM algorithm exclusively. Additionally, computation of the log-likelihood can be sped up by using parallel estimation via mirtCluster.

Usage

mixedmirt(
  data,
  covdata = NULL,
  model = 1,
  fixed = ~1,
  random = NULL,
  itemtype = "Rasch",
  lr.fixed = ~1,
  lr.random = NULL,
  itemdesign = NULL,
  constrain = NULL,
  pars = NULL,
  return.design = FALSE,
  SE = TRUE,
  internal_constraints = TRUE,
  technical = list(SEtol = 1e-04),
  ...
)

Arguments

data: a matrix or data.frame that consists of numerically ordered data, organized in the form of integers, with missing data coded as NA
covdata: a data.frame that consists of the nrow(data) by K 'person level' fixed and random predictors. If missing data are present in this object then the observations from covdata and data will be removed row-wise via the rowSums(is.na(covdata)) > 0
model: an object returned from, or a string to be passed to, mirt.model() to declare how the IRT model is to be estimated. See mirt.model and mirt for more detail
fixed: a right sided R formula for specifying the fixed effect (aka 'explanatory') predictors from covdata and itemdesign. To estimate the intercepts for each item the keyword items is reserved and automatically added to the itemdesign input. If any polytomous items are being model the items are argument is not valid since all intercept parameters are freely estimated and identified with the parameterizations found in mirt, and the first column in the fixed design matrix (commonly the intercept or a reference group) is omitted
random: a right sided formula or list of formulas containing crossed random effects of the form v1 + ... v_n | G, where G is the grouping variable and v_n are random numeric predictors within each group. If no intercept value is specified then by default the correlations between the v's and G are estimated, but can be suppressed by including the ~ -1 + ... or 0 constant. G may contain interaction terms, such as group:items to include cross or person-level interactions effects
itemtype: same as itemtype in mirt, except when the fixed or random inputs are used does not support the following item types: c('PC2PL', 'PC3PL', '2PLNRM', '3PLNRM', '3PLuNRM', '4PLNRM')
lr.fixed: an R formula (or list of formulas) to specify regression effects in the latent variables from the variables in covdata. This is used to construct models such as the so-called 'latent regression model' to explain person-level ability/trait differences. If a named list of formulas is supplied (where the names correspond to the latent trait names in model) then specific regression effects can be estimated for each factor. Supplying a single formula will estimate the regression parameters for all latent traits by default.
lr.random: a list of random effect terms for modeling variability in the latent trait scores, where the syntax uses the same style as in the random argument. Useful for building so-called 'multilevel IRT' models which are non-Rasch (multilevel Rasch models do not technically require these because they can be built using the fixed and random inputs alone)
itemdesign: a data.frame object used to create a design matrix for the items, where each nrow(itemdesign) == nitems and the number of columns is equal to the number of fixed effect predictors (i.e., item intercepts). By default an items variable is reserved for modeling the item intercept parameters
constrain: a list indicating parameter equality constrains. See mirt for more detail
pars: used for parameter starting values. See mirt for more detail
return.design: logical; return the design matrices before they have (potentially) been reassigned?
SE: logical; compute the standard errors by approximating the information matrix using the MHRM algorithm? Default is TRUE
internal_constraints: logical; use the internally defined constraints for constraining effects across persons and items? Default is TRUE. Setting this to FALSE runs the risk of under-identification
technical: the technical list passed to the MH-RM estimation engine, with the SEtol default increased to .0001. Additionally, the argument RANDSTART is available to indicate at which iteration (during the burn-in stage) the additional random effect variables should begin to be approximated (i.e., elements in lr.random and random). The default for RANDSTART is to start at iteration 100, and when random effects are included the default number of burn-in iterations is increased from 150 to 200. See mirt for further details
...: additional arguments to be passed to the MH-RM estimation engine. See mirt for more details and examples

Value

function returns an object of class MixedClass (MixedClass-class).

Details

For dichotomous response models, mixedmirt follows the general form

$$P(x = 1|\theta, \psi) = g + \frac{(u - g)}{1 + exp(-1 * [\theta a + X \beta + Z \delta])}$$

where X is a design matrix with associated $\beta$ fixed effect intercept coefficients, and Z is a design matrix with associated $\delta$ random effects for the intercepts. For simplicity and easier interpretation, the unique item intercept values typically found in $X \beta$ are extracted and reassigned within mirt's 'intercept' parameters (e.g., 'd'). To observe how the design matrices are structured prior to reassignment and estimation pass the argument return.design = TRUE.

Polytomous IRT models follow a similar format except the item intercepts are automatically estimated internally, rendering the items argument in the fixed formula redundant and therefore must be omitted from the specification. If there are a mixture of dichotomous and polytomous items the intercepts for the dichotomous models are also estimated for consistency.

The decomposition of the $\theta$ parameters is also possible to form latent regression and multilevel IRT models by using the lr.fixed and lr.random inputs. These effects decompose $\theta$ such that

$$\theta = V \Gamma + W \zeta + \epsilon$$

where V and W are fixed and random effects design matrices for the associated coefficients.

To simulate expected a posteriori predictions for the random effect terms use the randef function.

References

Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory Package for the R Environment. Journal of Statistical Software, 48(6), 1-29. doi:10.18637/jss.v048.i06

Chalmers, R. P. (2015). Extended Mixed-Effects Item Response Models with the MH-RM Algorithm. Journal of Educational Measurement, 52, 200-222. doi:10.1111/jedm.12072

Author

Phil Chalmers rphilip.chalmers@gmail.com

Examples


if (FALSE) { # \dontrun{

# make some data
set.seed(1234)
N <- 750
a <- matrix(rlnorm(10,.3,1),10,1)
d <- matrix(rnorm(10), 10)
Theta <- matrix(sort(rnorm(N)))
pseudoIQ <- Theta * 5 + 100  + rnorm(N, 0 , 5)
pseudoIQ <- (pseudoIQ - mean(pseudoIQ))/10  #rescale variable for numerical stability
group <- factor(rep(c('G1','G2','G3'), each = N/3))
data <- simdata(a,d,N, itemtype = rep('2PL',10), Theta=Theta)
covdata <- data.frame(group, pseudoIQ)

itemstats(data)

# use parallel computing
if(interactive()) mirtCluster()

# specify IRT model
model <- 'Theta = 1-10'

# model with no person predictors
mod0 <- mirt(data, model, itemtype = 'Rasch')

# group as a fixed effect predictor (aka, uniform dif)
mod1 <- mixedmirt(data, covdata, model, fixed = ~ 0 + group + items)
anova(mod0, mod1)
summary(mod1)
coef(mod1)

# same model as above in lme4
wide <- data.frame(id=1:nrow(data),data,covdata)
long <- reshape2::melt(wide, id.vars = c('id', 'group', 'pseudoIQ'))
library(lme4)
lmod0 <- glmer(value ~ 0 + variable + (1|id), long, family = binomial)
lmod1 <- glmer(value ~ 0 + group + variable + (1|id), long, family = binomial)
anova(lmod0, lmod1)

# model using 2PL items instead of Rasch
mod1b <- mixedmirt(data, covdata, model, fixed = ~ 0 + group + items, itemtype = '2PL')
anova(mod1, mod1b) #better with 2PL models using all criteria (as expected, given simdata pars)

# continuous predictor with group
mod2 <- mixedmirt(data, covdata, model, fixed = ~ 0 + group + items + pseudoIQ)
summary(mod2)
anova(mod1b, mod2)

# view fixed design matrix with and without unique item level intercepts
withint <- mixedmirt(data, covdata, model, fixed = ~ 0 + items + group, return.design = TRUE)
withoutint <- mixedmirt(data, covdata, model, fixed = ~ 0 + group, return.design = TRUE)

# notice that in result above, the intercepts 'items1 to items 10' were reassigned to 'd'
head(withint$X)
tail(withint$X)
head(withoutint$X) # no intercepts design here to be reassigned into item intercepts
tail(withoutint$X)

###################################################
### random effects
# make the number of groups much larger
covdata$group <- factor(rep(paste0('G',1:50), each = N/50))

# random groups
rmod1 <- mixedmirt(data, covdata, 1, fixed = ~ 0 + items, random = ~ 1|group)
summary(rmod1)
coef(rmod1)

# random groups and random items
rmod2 <- mixedmirt(data, covdata, 1, random = list(~ 1|group, ~ 1|items))
summary(rmod2)
eff <- randef(rmod2) #estimate random effects

# random slopes with fixed intercepts (suppressed correlation)
rmod3 <- mixedmirt(data, covdata, 1, fixed = ~ 0 + items, random = ~ -1 + pseudoIQ|group)
summary(rmod3)
eff <- randef(rmod3)
str(eff)

###################################################
## LLTM, and 2PL version of LLTM
data(SAT12)
data <- key2binary(SAT12,
                   key = c(1,4,5,2,3,1,2,1,3,1,2,4,2,1,5,3,4,4,1,4,3,3,4,1,3,5,1,3,1,5,4,5))
model <- 'Theta = 1-32'

# for unconditional intercept comparison
mod <- mirt(data, model, itemtype='Rasch')
coef(mod, simplify=TRUE)

# Suppose that the first 16 items were suspected to be easier than the last 16 items,
#   and we wish to test this item structure hypothesis (more intercept designs are possible
#   by including more columns).
itemdesign <- data.frame(itemorder = factor(c(rep('easier', 16), rep('harder', 16))))
rownames(itemdesign) <- colnames(data)
itemdesign

# notice that the 'fixed = ~ ... + items' argument is omitted
LLTM <- mixedmirt(data, model = model, fixed = ~ 0 + itemorder, itemdesign = itemdesign,
   SE = TRUE) # SE argument ensures that the information matrix is computed accurately
summary(LLTM)
coef(LLTM)
wald(LLTM)
L <- matrix(c(-1, 1, 0), 1)
wald(LLTM, L) #first half different from second

# compare to items with estimated slopes (2PL)
twoPL <- mixedmirt(data, model = model, fixed = ~ 0 + itemorder, itemtype = '2PL',
                   itemdesign = itemdesign)
# twoPL not mixing too well (AR should be between .2 and .5), decrease MHcand
twoPL <- mixedmirt(data, model = model, fixed = ~ 0 + itemorder, itemtype = '2PL',
                  itemdesign = itemdesign, technical = list(MHcand = 0.8))
anova(twoPL, LLTM) #much better fit
summary(twoPL)
coef(twoPL)

wald(twoPL)
L <- matrix(0, 1, 34)
L[1, 1] <- 1
L[1, 2] <- -1
wald(twoPL, L) # n.s., which is the correct conclusion. Rasch approach gave wrong inference

## LLTM with item error term
LLTMwithError <- mixedmirt(data, model = model, fixed = ~ 0 + itemorder, random = ~ 1|items,
    itemdesign = itemdesign)
summary(LLTMwithError)
# large item level variance after itemorder is regressed; not a great predictor of item difficulty
coef(LLTMwithError)

###################################################
### Polytomous example

# make an arbitrary group difference
covdat <- data.frame(group = rep(c('m', 'f'), nrow(Science)/2))

# partial credit model
mod <- mixedmirt(Science, covdat, model=1, fixed = ~ 0 + group)
coef(mod)

# gpcm to estimate slopes
mod2 <- mixedmirt(Science, covdat, model=1, fixed = ~ 0 + group,
                 itemtype = 'gpcm')
summary(mod2)
anova(mod, mod2)

# graded model
mod3 <- mixedmirt(Science, covdat, model=1, fixed = ~ 0 + group,
                 itemtype = 'graded')
coef(mod3)


###################################################
# latent regression with Rasch and 2PL models

set.seed(1)
n <- 300
a <- matrix(1, 10)
d <- matrix(rnorm(10))
Theta <- matrix(c(rnorm(n, 0), rnorm(n, 1), rnorm(n, 2)))
covdata <- data.frame(group=rep(c('g1','g2','g3'), each=n))
dat <- simdata(a, d, N=n*3, Theta=Theta, itemtype = '2PL')
itemstats(dat)

# had we known the latent abilities, we could have computed the regression coefs
summary(lm(Theta ~ covdata$group))

# but all we have is observed test data. Latent regression helps to recover these coefs
# Rasch model approach (and mirt equivalent)
rmod0 <- mirt(dat, 1, 'Rasch') # unconditional

# these two models are equivalent
rmod1a <- mirt(dat, 1, 'Rasch', covdata = covdata, formula = ~ group)
rmod1b <- mixedmirt(dat, covdata, 1, fixed = ~ 0 + items + group)
anova(rmod0, rmod1b)
coef(rmod1a, simplify=TRUE)
summary(rmod1b)

# 2PL, requires different input to allow Theta variance to remain fixed
mod0 <- mirt(dat, 1) # unconditional
mod1a <- mirt(dat, 1, covdata = covdata, formula = ~ group, itemtype = '2PL')
mod1b <- mixedmirt(dat, covdata, 1, fixed = ~ 0 + items, lr.fixed = ~group, itemtype = '2PL')
anova(mod0, mod1b)
coef(mod1a)$lr.betas
summary(mod1b)

# specifying specific regression effects is accomplished by passing a list of formula
model <- 'F1 = 1-5
         F2 = 6-10'
covdata$contvar <- rnorm(nrow(covdata))
mod2 <- mirt(dat, model, itemtype = 'Rasch', covdata=covdata,
        formula = list(F1 = ~ group + contvar, F2 = ~ group))
coef(mod2)[11:12]
mod2b <- mixedmirt(dat, covdata, model, fixed = ~ 0 + items,
        lr.fixed = list(F1 = ~ group + contvar, F2 = ~ group))
summary(mod2b)

####################################################
## Simulated Multilevel Rasch Model

set.seed(1)
N <- 2000
a <- matrix(rep(1,10),10,1)
d <- matrix(rnorm(10))
cluster = 100
random_intercept = rnorm(cluster,0,1)
Theta = numeric()
for (i in 1:cluster)
    Theta <- c(Theta, rnorm(N/cluster,0,1) + random_intercept[i])

group = factor(rep(paste0('G',1:cluster), each = N/cluster))
covdata <- data.frame(group)
dat <- simdata(a,d,N, itemtype = rep('2PL',10), Theta=matrix(Theta))
itemstats(dat)

# null model
mod1 <- mixedmirt(dat, covdata, 1, fixed = ~ 0 + items, random = ~ 1|group)
summary(mod1)

# include level 2 predictor for 'group' variance
covdata$group_pred <- rep(random_intercept, each = N/cluster)
mod2 <- mixedmirt(dat, covdata, 1, fixed = ~ 0 + items + group_pred, random = ~ 1|group)

# including group means predicts nearly all variability in 'group'
summary(mod2)
anova(mod1, mod2)

# can also be fit for Rasch/non-Rasch models with the lr.random input
mod1b <- mixedmirt(dat, covdata, 1, fixed = ~ 0 + items, lr.random = ~ 1|group)
summary(mod1b)

mod2b <- mixedmirt(dat, covdata, 1, fixed = ~ 0 + items + group_pred, lr.random = ~ 1|group)
summary(mod2b)
anova(mod1b, mod2b)

mod3 <- mixedmirt(dat, covdata, 1, fixed = ~ 0 + items, lr.random = ~ 1|group, itemtype = '2PL')
summary(mod3)
anova(mod1b, mod3)

head(cbind(randef(mod3)$group, random_intercept))

} # }