R: Mixed effects modeling for MIRT models

mixedmirt {mirt}

R Documentation

Mixed effects modeling for MIRT models

Description

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, 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

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.

Value

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

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

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

Examples

## No test: 

# 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)

## $overall
##    N mean_total.score sd_total.score ave.r  sd.r alpha
##  750            4.655          2.346 0.166 0.133 0.671
## 
## $itemstats
##           N  mean    sd total.r total.r_if_rm alpha_if_rm
## Item_1  750 0.363 0.481   0.368         0.172       0.678
## Item_2  750 0.335 0.472   0.631         0.485       0.617
## Item_3  750 0.428 0.495   0.711         0.580       0.594
## Item_4  750 0.512 0.500   0.234         0.021       0.708
## Item_5  750 0.628 0.484   0.638         0.490       0.615
## Item_6  750 0.472 0.500   0.669         0.523       0.607
## Item_7  750 0.385 0.487   0.471         0.286       0.657
## Item_8  750 0.301 0.459   0.483         0.312       0.651
## Item_9  750 0.319 0.466   0.481         0.307       0.652
## Item_10 750 0.912 0.283   0.295         0.180       0.670
## 
## $proportions
##             0     1
## Item_1  0.637 0.363
## Item_2  0.665 0.335
## Item_3  0.572 0.428
## Item_4  0.488 0.512
## Item_5  0.372 0.628
## Item_6  0.528 0.472
## Item_7  0.615 0.385
## Item_8  0.699 0.301
## Item_9  0.681 0.319
## Item_10 0.088 0.912

# 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)

##           AIC    SABIC       HQ      BIC    logLik      X2 df p
## mod0 8799.543 8815.435 8819.126 8850.364 -4388.772             
## mod1 8111.225 8130.006 8134.368 8171.286 -4042.613 692.318  2 0

summary(mod1)

## 
## Call:
## mixedmirt(data = data, covdata = covdata, model = model, fixed = ~0 + 
##     group + items)
## 
## --------------
## FIXED EFFECTS:
##         Estimate Std.Error z.value
## groupG1   -1.858     0.100 -18.495
## groupG2   -0.748     0.094  -7.979
## groupG3    0.515     0.094   5.493
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##       Theta
## Theta 0.118

coef(mod1)

## $Item_1
##         groupG1 groupG2 groupG3 a1  d  g  u
## par      -1.858  -0.748   0.515  1  0  0  1
## CI_2.5   -2.055  -0.931   0.331 NA NA NA NA
## CI_97.5  -1.661  -0.564   0.699 NA NA NA NA
## 
## $Item_2
##         groupG1 groupG2 groupG3 a1      d  g  u
## par      -1.858  -0.748   0.515  1 -0.152  0  1
## CI_2.5   -2.055  -0.931   0.331 NA -0.388 NA NA
## CI_97.5  -1.661  -0.564   0.699 NA  0.084 NA NA
## 
## $Item_3
##         groupG1 groupG2 groupG3 a1     d  g  u
## par      -1.858  -0.748   0.515  1 0.340  0  1
## CI_2.5   -2.055  -0.931   0.331 NA 0.109 NA NA
## CI_97.5  -1.661  -0.564   0.699 NA 0.572 NA NA
## 
## $Item_4
##         groupG1 groupG2 groupG3 a1     d  g  u
## par      -1.858  -0.748   0.515  1 0.763  0  1
## CI_2.5   -2.055  -0.931   0.331 NA 0.532 NA NA
## CI_97.5  -1.661  -0.564   0.699 NA 0.994 NA NA
## 
## $Item_5
##         groupG1 groupG2 groupG3 a1     d  g  u
## par      -1.858  -0.748   0.515  1 1.353  0  1
## CI_2.5   -2.055  -0.931   0.331 NA 1.119 NA NA
## CI_97.5  -1.661  -0.564   0.699 NA 1.588 NA NA
## 
## $Item_6
##         groupG1 groupG2 groupG3 a1     d  g  u
## par      -1.858  -0.748   0.515  1 0.563  0  1
## CI_2.5   -2.055  -0.931   0.331 NA 0.332 NA NA
## CI_97.5  -1.661  -0.564   0.699 NA 0.793 NA NA
## 
## $Item_7
##         groupG1 groupG2 groupG3 a1      d  g  u
## par      -1.858  -0.748   0.515  1  0.120  0  1
## CI_2.5   -2.055  -0.931   0.331 NA -0.113 NA NA
## CI_97.5  -1.661  -0.564   0.699 NA  0.353 NA NA
## 
## $Item_8
##         groupG1 groupG2 groupG3 a1      d  g  u
## par      -1.858  -0.748   0.515  1 -0.339  0  1
## CI_2.5   -2.055  -0.931   0.331 NA -0.578 NA NA
## CI_97.5  -1.661  -0.564   0.699 NA -0.101 NA NA
## 
## $Item_9
##         groupG1 groupG2 groupG3 a1      d  g  u
## par      -1.858  -0.748   0.515  1 -0.241  0  1
## CI_2.5   -2.055  -0.931   0.331 NA -0.478 NA NA
## CI_97.5  -1.661  -0.564   0.699 NA -0.004 NA NA
## 
## $Item_10
##         groupG1 groupG2 groupG3 a1     d  g  u
## par      -1.858  -0.748   0.515  1 3.432  0  1
## CI_2.5   -2.055  -0.931   0.331 NA 3.117 NA NA
## CI_97.5  -1.661  -0.564   0.699 NA 3.747 NA NA
## 
## $GroupPars
##         MEAN_1 COV_11
## par          0  0.118
## CI_2.5      NA  0.072
## CI_97.5     NA  0.164

# 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)

## Loading required package: Matrix

## 
## Attaching package: 'lme4'

## The following object is masked from 'package:mirt':
## 
##     fixef

lmod0 <- glmer(value ~ 0 + variable + (1|id), long, family = binomial)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0292334 (tol = 0.002, component 1)

lmod1 <- glmer(value ~ 0 + group + variable + (1|id), long, family = binomial)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00585234 (tol = 0.002, component 1)

anova(lmod0, lmod1)

## Data: long
## Models:
## lmod0: value ~ 0 + variable + (1 | id)
## lmod1: value ~ 0 + group + variable + (1 | id)
##       npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## lmod0   11 8813.5 8889.6 -4395.7   8791.5                         
## lmod1   13 8109.5 8199.5 -4041.8   8083.5 707.93  2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 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)

##            AIC    SABIC       HQ      BIC    logLik      X2 df p
## mod1  8111.225 8130.006 8134.368 8171.286 -4042.613             
## mod1b 7973.762 8005.545 8012.926 8075.403 -3964.881 155.463  9 0

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

## 
## Call:
## mixedmirt(data = data, covdata = covdata, model = model, fixed = ~0 + 
##     group + items + pseudoIQ)
## 
## --------------
## FIXED EFFECTS:
##          Estimate Std.Error z.value
## groupG1    -1.711     0.105 -16.317
## groupG2    -0.750     0.093  -8.029
## groupG3     0.366     0.099   3.714
## pseudoIQ    0.268     0.058   4.594
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##       Theta
## Theta 0.105

anova(mod1b, mod2)

##            AIC    SABIC       HQ      BIC    logLik       X2 df   p
## mod1b 7973.762 8005.545 8012.926 8075.403 -3964.881                
## mod2  8089.308 8109.533 8114.231 8153.989 -4030.654 -131.546 -8 NaN

# 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)

##     items1 items2 items3 items4 items5 items6 items7 items8 items9 items10
## 1.1      1      0      0      0      0      0      0      0      0       0
## 2.1      1      0      0      0      0      0      0      0      0       0
## 3.1      1      0      0      0      0      0      0      0      0       0
## 4.1      1      0      0      0      0      0      0      0      0       0
## 5.1      1      0      0      0      0      0      0      0      0       0
## 6.1      1      0      0      0      0      0      0      0      0       0
##     groupG2 groupG3
## 1.1       0       0
## 2.1       0       0
## 3.1       0       0
## 4.1       0       0
## 5.1       0       0
## 6.1       0       0

tail(withint$X)

##        items1 items2 items3 items4 items5 items6 items7 items8 items9 items10
## 745.10      0      0      0      0      0      0      0      0      0       1
## 746.10      0      0      0      0      0      0      0      0      0       1
## 747.10      0      0      0      0      0      0      0      0      0       1
## 748.10      0      0      0      0      0      0      0      0      0       1
## 749.10      0      0      0      0      0      0      0      0      0       1
## 750.10      0      0      0      0      0      0      0      0      0       1
##        groupG2 groupG3
## 745.10       0       1
## 746.10       0       1
## 747.10       0       1
## 748.10       0       1
## 749.10       0       1
## 750.10       0       1

head(withoutint$X) # no intercepts design here to be reassigned into item intercepts

##     groupG1 groupG2 groupG3
## 1.1       1       0       0
## 2.1       1       0       0
## 3.1       1       0       0
## 4.1       1       0       0
## 5.1       1       0       0
## 6.1       1       0       0

tail(withoutint$X)

##        groupG1 groupG2 groupG3
## 745.10       0       0       1
## 746.10       0       0       1
## 747.10       0       0       1
## 748.10       0       0       1
## 749.10       0       0       1
## 750.10       0       0       1

###################################################
### 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)

## 
## Call:
## mixedmirt(data = data, covdata = covdata, model = 1, fixed = ~0 + 
##     items, random = ~1 | group)
## 
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##       F1
## F1 0.062
## 
## $group
##           COV_group
## COV_group      1.11

coef(rmod1)

## $Item_1
##         a1      d  g  u
## par      1 -0.713  0  1
## CI_2.5  NA -0.927 NA NA
## CI_97.5 NA -0.500 NA NA
## 
## $Item_2
##         a1      d  g  u
## par      1 -0.867  0  1
## CI_2.5  NA -1.081 NA NA
## CI_97.5 NA -0.653 NA NA
## 
## $Item_3
##         a1      d  g  u
## par      1 -0.370  0  1
## CI_2.5  NA -0.583 NA NA
## CI_97.5 NA -0.156 NA NA
## 
## $Item_4
##         a1      d  g  u
## par      1  0.057  0  1
## CI_2.5  NA -0.161 NA NA
## CI_97.5 NA  0.276 NA NA
## 
## $Item_5
##         a1     d  g  u
## par      1 0.654  0  1
## CI_2.5  NA 0.422 NA NA
## CI_97.5 NA 0.886 NA NA
## 
## $Item_6
##         a1      d  g  u
## par      1 -0.145  0  1
## CI_2.5  NA -0.361 NA NA
## CI_97.5 NA  0.071 NA NA
## 
## $Item_7
##         a1      d  g  u
## par      1 -0.592  0  1
## CI_2.5  NA -0.805 NA NA
## CI_97.5 NA -0.379 NA NA
## 
## $Item_8
##         a1      d  g  u
## par      1 -1.057  0  1
## CI_2.5  NA -1.273 NA NA
## CI_97.5 NA -0.842 NA NA
## 
## $Item_9
##         a1      d  g  u
## par      1 -0.957  0  1
## CI_2.5  NA -1.172 NA NA
## CI_97.5 NA -0.743 NA NA
## 
## $Item_10
##         a1     d  g  u
## par      1 2.769  0  1
## CI_2.5  NA 2.425 NA NA
## CI_97.5 NA 3.113 NA NA
## 
## $GroupPars
##         MEAN_1 COV_11
## par          0  0.062
## CI_2.5      NA  0.004
## CI_97.5     NA  0.120
## 
## $group
##         COV_group_group
## par               1.112
## CI_2.5            0.631
## CI_97.5           1.592

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

## 
## Call:
## mixedmirt(data = data, covdata = covdata, model = 1, random = list(~1 | 
##     group, ~1 | items))
## 
## --------------
## FIXED EFFECTS:
##             Estimate Std.Error z.value
## (Intercept)   -0.605     0.013 -47.455
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##        F1
## F1 0.0593
## 
## $group
##           COV_group
## COV_group     0.969
## 
## $items
##           COV_items
## COV_items     0.694

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)

## 
## Call:
## mixedmirt(data = data, covdata = covdata, model = 1, fixed = ~0 + 
##     items, random = ~-1 + pseudoIQ | group)
## 
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##       F1
## F1 0.073
## 
## $group
##              COV_group COV_pseudoIQ
## COV_group        0.903         0.00
## COV_pseudoIQ     0.000         0.16

eff <- randef(rmod3)
str(eff)

## List of 2
##  $ Theta: num [1:750, 1] -0.0706 0.0482 0.0494 0.0601 -0.0411 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "F1"
##  $ group: num [1:50, 1:2] -1.68 -1.18 -1.33 -1.6 -1.41 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:50] "G1" "G2" "G3" "G4" ...
##   .. ..$ : chr [1:2] "group" "pseudoIQ"

###################################################
## 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'

# 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))))

# 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)

## 
## Call:
## mixedmirt(data = data, model = model, fixed = ~0 + itemorder, 
##     itemdesign = itemdesign, SE = TRUE)
## 
## --------------
## FIXED EFFECTS:
##                 Estimate Std.Error z.value
## itemordereasier    0.165     0.029   5.746
## itemorderharder    0.456     0.029  15.757
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##       Theta
## Theta 0.359

coef(LLTM)

## $Item.1
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.2
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.3
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.4
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.5
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.6
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.7
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.8
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.9
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.10
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.11
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.12
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.13
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.14
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.15
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.16
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.17
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.18
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.19
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.20
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.21
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.22
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.23
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.24
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.25
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.26
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.27
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.28
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.29
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.30
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.31
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $Item.32
##         itemordereasier itemorderharder a1  d  g  u
## par               0.165           0.456  1  0  0  1
## CI_2.5            0.109           0.400 NA NA NA NA
## CI_97.5           0.221           0.513 NA NA NA NA
## 
## $GroupPars
##         MEAN_1 COV_11
## par          0  0.359
## CI_2.5      NA  0.300
## CI_97.5     NA  0.417

wald(LLTM)

## itemordereasier.1.7.13.19.25.31.37.43.49.55.61.67.73.79.85.91.97.103.109.115.121.127.133.139.145.151.157.163.169.175.181.187 
##                                                                                                                        0.165 
## itemorderharder.2.8.14.20.26.32.38.44.50.56.62.68.74.80.86.92.98.104.110.116.122.128.134.140.146.152.158.164.170.176.182.188 
##                                                                                                                        0.456 
##                                                                                                                   COV_11.194 
##                                                                                                                        0.359

L <- matrix(c(-1, 1, 0), 1)
wald(LLTM, L) #first half different from second

##          W df p
## 1 92.08467  1 0

# 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

##            AIC    SABIC       HQ      BIC   logLik       X2  df   p
## twoPL 20479.20 20520.76 20537.40 20628.70 -10205.6                 
## LLTM  25463.99 25467.66 25469.13 25477.19 -12729.0 -5046.79 -31 NaN

summary(twoPL)

## 
## Call:
## mixedmirt(data = data, model = model, fixed = ~0 + itemorder, 
##     itemtype = "2PL", itemdesign = itemdesign, technical = list(MHcand = 0.8))
## 
## --------------
## FIXED EFFECTS:
##                 Estimate Std.Error z.value
## itemordereasier   -1.669     0.087 -19.127
## itemorderharder   -1.644     0.095 -17.316
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##       Theta
## Theta     1

coef(twoPL)

## $Item.1
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 0.930  0  0  1
## CI_2.5           -1.840          -1.831 0.696 NA NA NA
## CI_97.5          -1.498          -1.458 1.165 NA NA NA
## 
## $Item.2
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 2.526  0  0  1
## CI_2.5           -1.840          -1.831 2.201 NA NA NA
## CI_97.5          -1.498          -1.458 2.851 NA NA NA
## 
## $Item.3
##         itemordereasier itemorderharder   a1  d  g  u
## par              -1.669          -1.644 1.00  0  0  1
## CI_2.5           -1.840          -1.831 0.77 NA NA NA
## CI_97.5          -1.498          -1.458 1.23 NA NA NA
## 
## $Item.4
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 1.288  0  0  1
## CI_2.5           -1.840          -1.831 1.035 NA NA NA
## CI_97.5          -1.498          -1.458 1.540 NA NA NA
## 
## $Item.5
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 2.704  0  0  1
## CI_2.5           -1.840          -1.831 2.359 NA NA NA
## CI_97.5          -1.498          -1.458 3.049 NA NA NA
## 
## $Item.6
##         itemordereasier itemorderharder   a1  d  g  u
## par              -1.669          -1.644 0.42  0  0  1
## CI_2.5           -1.840          -1.831 0.19 NA NA NA
## CI_97.5          -1.498          -1.458 0.65 NA NA NA
## 
## $Item.7
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 3.913  0  0  1
## CI_2.5           -1.840          -1.831 3.461 NA NA NA
## CI_97.5          -1.498          -1.458 4.365 NA NA NA
## 
## $Item.8
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 0.513  0  0  1
## CI_2.5           -1.840          -1.831 0.282 NA NA NA
## CI_97.5          -1.498          -1.458 0.744 NA NA NA
## 
## $Item.9
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 5.455  0  0  1
## CI_2.5           -1.840          -1.831 4.840 NA NA NA
## CI_97.5          -1.498          -1.458 6.071 NA NA NA
## 
## $Item.10
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 1.628  0  0  1
## CI_2.5           -1.840          -1.831 1.352 NA NA NA
## CI_97.5          -1.498          -1.458 1.905 NA NA NA
## 
## $Item.11
##         itemordereasier itemorderharder     a1  d  g  u
## par              -1.669          -1.644 13.512  0  0  1
## CI_2.5           -1.840          -1.831 10.680 NA NA NA
## CI_97.5          -1.498          -1.458 16.344 NA NA NA
## 
## $Item.12
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 1.245  0  0  1
## CI_2.5           -1.840          -1.831 0.974 NA NA NA
## CI_97.5          -1.498          -1.458 1.515 NA NA NA
## 
## $Item.13
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 3.042  0  0  1
## CI_2.5           -1.840          -1.831 2.668 NA NA NA
## CI_97.5          -1.498          -1.458 3.416 NA NA NA
## 
## $Item.14
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 3.488  0  0  1
## CI_2.5           -1.840          -1.831 3.076 NA NA NA
## CI_97.5          -1.498          -1.458 3.899 NA NA NA
## 
## $Item.15
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 4.749  0  0  1
## CI_2.5           -1.840          -1.831 4.225 NA NA NA
## CI_97.5          -1.498          -1.458 5.272 NA NA NA
## 
## $Item.16
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 1.469  0  0  1
## CI_2.5           -1.840          -1.831 1.205 NA NA NA
## CI_97.5          -1.498          -1.458 1.733 NA NA NA
## 
## $Item.17
##         itemordereasier itemorderharder     a1  d  g  u
## par              -1.669          -1.644  9.952  0  0  1
## CI_2.5           -1.840          -1.831  8.451 NA NA NA
## CI_97.5          -1.498          -1.458 11.454 NA NA NA
## 
## $Item.18
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 1.408  0  0  1
## CI_2.5           -1.840          -1.831 1.145 NA NA NA
## CI_97.5          -1.498          -1.458 1.670 NA NA NA
## 
## $Item.19
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 2.155  0  0  1
## CI_2.5           -1.840          -1.831 1.825 NA NA NA
## CI_97.5          -1.498          -1.458 2.486 NA NA NA
## 
## $Item.20
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 5.787  0  0  1
## CI_2.5           -1.840          -1.831 5.118 NA NA NA
## CI_97.5          -1.498          -1.458 6.456 NA NA NA
## 
## $Item.21
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 6.345  0  0  1
## CI_2.5           -1.840          -1.831 5.595 NA NA NA
## CI_97.5          -1.498          -1.458 7.094 NA NA NA
## 
## $Item.22
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 7.826  0  0  1
## CI_2.5           -1.840          -1.831 6.827 NA NA NA
## CI_97.5          -1.498          -1.458 8.824 NA NA NA
## 
## $Item.23
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 0.982  0  0  1
## CI_2.5           -1.840          -1.831 0.732 NA NA NA
## CI_97.5          -1.498          -1.458 1.231 NA NA NA
## 
## $Item.24
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 3.622  0  0  1
## CI_2.5           -1.840          -1.831 3.186 NA NA NA
## CI_97.5          -1.498          -1.458 4.057 NA NA NA
## 
## $Item.25
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 1.293  0  0  1
## CI_2.5           -1.840          -1.831 1.021 NA NA NA
## CI_97.5          -1.498          -1.458 1.565 NA NA NA
## 
## $Item.26
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 1.907  0  0  1
## CI_2.5           -1.840          -1.831 1.608 NA NA NA
## CI_97.5          -1.498          -1.458 2.205 NA NA NA
## 
## $Item.27
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 5.656  0  0  1
## CI_2.5           -1.840          -1.831 5.018 NA NA NA
## CI_97.5          -1.498          -1.458 6.293 NA NA NA
## 
## $Item.28
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 2.120  0  0  1
## CI_2.5           -1.840          -1.831 1.807 NA NA NA
## CI_97.5          -1.498          -1.458 2.433 NA NA NA
## 
## $Item.29
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 1.165  0  0  1
## CI_2.5           -1.840          -1.831 0.910 NA NA NA
## CI_97.5          -1.498          -1.458 1.420 NA NA NA
## 
## $Item.30
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 1.442  0  0  1
## CI_2.5           -1.840          -1.831 1.156 NA NA NA
## CI_97.5          -1.498          -1.458 1.727 NA NA NA
## 
## $Item.31
##         itemordereasier itemorderharder    a1  d  g  u
## par              -1.669          -1.644 5.235  0  0  1
## CI_2.5           -1.840          -1.831 4.637 NA NA NA
## CI_97.5          -1.498          -1.458 5.832 NA NA NA
## 
## $Item.32
##         itemordereasier itemorderharder     a1  d  g  u
## par              -1.669          -1.644  0.097  0  0  1
## CI_2.5           -1.840          -1.831 -0.163 NA NA NA
## CI_97.5          -1.498          -1.458  0.357 NA NA NA
## 
## $GroupPars
##         MEAN_1 COV_11
## par          0      1
## CI_2.5      NA     NA
## CI_97.5     NA     NA

wald(twoPL)

## itemordereasier.1.7.13.19.25.31.37.43.49.55.61.67.73.79.85.91.97.103.109.115.121.127.133.139.145.151.157.163.169.175.181.187 
##                                                                                                                       -1.669 
## itemorderharder.2.8.14.20.26.32.38.44.50.56.62.68.74.80.86.92.98.104.110.116.122.128.134.140.146.152.158.164.170.176.182.188 
##                                                                                                                       -1.644 
##                                                                                                                         a1.3 
##                                                                                                                        0.930 
##                                                                                                                         a1.9 
##                                                                                                                        2.526 
##                                                                                                                        a1.15 
##                                                                                                                        1.000 
##                                                                                                                        a1.21 
##                                                                                                                        1.288 
##                                                                                                                        a1.27 
##                                                                                                                        2.704 
##                                                                                                                        a1.33 
##                                                                                                                        0.420 
##                                                                                                                        a1.39 
##                                                                                                                        3.913 
##                                                                                                                        a1.45 
##                                                                                                                        0.513 
##                                                                                                                        a1.51 
##                                                                                                                        5.455 
##                                                                                                                        a1.57 
##                                                                                                                        1.628 
##                                                                                                                        a1.63 
##                                                                                                                       13.512 
##                                                                                                                        a1.69 
##                                                                                                                        1.245 
##                                                                                                                        a1.75 
##                                                                                                                        3.042 
##                                                                                                                        a1.81 
##                                                                                                                        3.488 
##                                                                                                                        a1.87 
##                                                                                                                        4.749 
##                                                                                                                        a1.93 
##                                                                                                                        1.469 
##                                                                                                                        a1.99 
##                                                                                                                        9.952 
##                                                                                                                       a1.105 
##                                                                                                                        1.408 
##                                                                                                                       a1.111 
##                                                                                                                        2.155 
##                                                                                                                       a1.117 
##                                                                                                                        5.787 
##                                                                                                                       a1.123 
##                                                                                                                        6.345 
##                                                                                                                       a1.129 
##                                                                                                                        7.826 
##                                                                                                                       a1.135 
##                                                                                                                        0.982 
##                                                                                                                       a1.141 
##                                                                                                                        3.622 
##                                                                                                                       a1.147 
##                                                                                                                        1.293 
##                                                                                                                       a1.153 
##                                                                                                                        1.907 
##                                                                                                                       a1.159 
##                                                                                                                        5.656 
##                                                                                                                       a1.165 
##                                                                                                                        2.120 
##                                                                                                                       a1.171 
##                                                                                                                        1.165 
##                                                                                                                       a1.177 
##                                                                                                                        1.442 
##                                                                                                                       a1.183 
##                                                                                                                        5.235 
##                                                                                                                       a1.189 
##                                                                                                                        0.097

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

##            W df         p
## 1 0.07451603  1 0.7848715

## LLTM with item error term
LLTMwithError <- mixedmirt(data, model = model, fixed = ~ 0 + itemorder, random = ~ 1|items,
    itemdesign = itemdesign)
summary(LLTMwithError)

## 
## Call:
## mixedmirt(data = data, model = model, fixed = ~0 + itemorder, 
##     random = ~1 | items, itemdesign = itemdesign)
## 
## --------------
## FIXED EFFECTS:
##                 Estimate Std.Error z.value
## itemordereasier    0.147     0.153   0.961
## itemorderharder    0.567     0.138   4.107
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##       Theta
## Theta  0.77
## 
## $items
##           COV_items
## COV_items      2.36

# large item level variance after itemorder is regressed; not a great predictor of item difficulty
coef(LLTMwithError)

## $Item.1
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.2
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.3
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.4
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.5
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.6
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.7
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.8
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.9
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.10
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.11
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.12
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.13
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.14
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.15
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.16
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.17
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.18
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.19
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.20
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.21
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.22
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.23
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.24
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.25
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.26
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.27
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.28
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.29
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.30
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.31
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $Item.32
##         itemordereasier itemorderharder a1  d  g  u
## par               0.147           0.567  1  0  0  1
## CI_2.5           -0.153           0.296 NA NA NA NA
## CI_97.5           0.446           0.838 NA NA NA NA
## 
## $GroupPars
##         MEAN_1 COV_11
## par          0  0.770
## CI_2.5      NA  0.656
## CI_97.5     NA  0.884
## 
## $items
##         COV_items_items
## par               2.363
## CI_2.5            1.229
## CI_97.5           3.497

###################################################
### 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)

## $Comfort
##         groupm a1 ak0 ak1 ak2 ak3 d0    d1    d2    d3
## par     -0.084  1   0   1   2   3  0 3.098 5.718 4.364
## CI_2.5  -0.336 NA  NA  NA  NA  NA NA 2.086 4.677 3.257
## CI_97.5  0.167 NA  NA  NA  NA  NA NA 4.110 6.759 5.472
## 
## $Work
##         groupm a1 ak0 ak1 ak2 ak3 d0    d1    d2    d3
## par     -0.084  1   0   1   2   3  0 1.919 2.859 1.038
## CI_2.5  -0.336 NA  NA  NA  NA  NA NA 1.455 2.311 0.352
## CI_97.5  0.167 NA  NA  NA  NA  NA NA 2.382 3.406 1.723
## 
## $Future
##         groupm a1 ak0 ak1 ak2 ak3 d0    d1    d2    d3
## par     -0.084  1   0   1   2   3  0 2.665 4.114 3.013
## CI_2.5  -0.336 NA  NA  NA  NA  NA NA 2.026 3.406 2.214
## CI_97.5  0.167 NA  NA  NA  NA  NA NA 3.304 4.822 3.813
## 
## $Benefit
##         groupm a1 ak0 ak1 ak2 ak3 d0    d1    d2    d3
## par     -0.084  1   0   1   2   3  0 2.469 3.398 2.076
## CI_2.5  -0.336 NA  NA  NA  NA  NA NA 1.932 2.779 1.349
## CI_97.5  0.167 NA  NA  NA  NA  NA NA 3.006 4.016 2.803
## 
## $GroupPars
##         MEAN_1 COV_11
## par          0  0.985
## CI_2.5      NA  0.707
## CI_97.5     NA  1.264

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

## 
## Call:
## mixedmirt(data = Science, covdata = covdat, model = 1, fixed = ~0 + 
##     group, itemtype = "gpcm")
## 
## --------------
## FIXED EFFECTS:
##        Estimate Std.Error z.value
## groupm   -0.176     0.118  -1.494
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##    F1
## F1  1

anova(mod, mod2)

##           AIC    SABIC       HQ      BIC    logLik     X2 df     p
## mod  3268.671 3279.847 3290.706 3324.269 -1620.336                
## mod2 3258.272 3271.844 3285.029 3325.784 -1612.136 16.399  3 0.001

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

## $Comfort
##         groupm    a1    d1    d2     d3
## par     -0.248 1.000 4.951 2.735 -1.327
## CI_2.5  -0.578 0.636 3.989 2.289 -1.682
## CI_97.5  0.082 1.363 5.913 3.181 -0.971
## 
## $Work
##         groupm    a1    d1    d2     d3
## par     -0.248 1.212 3.034 1.020 -2.133
## CI_2.5  -0.578 0.873 2.545 0.697 -2.554
## CI_97.5  0.082 1.551 3.522 1.342 -1.712
## 
## $Future
##         groupm    a1    d1    d2     d3
## par     -0.248 2.581 5.761 2.524 -1.995
## CI_2.5  -0.578 1.251 3.733 1.497 -2.807
## CI_97.5  0.082 3.911 7.789 3.551 -1.182
## 
## $Benefit
##         groupm    a1    d1    d2     d3
## par     -0.248 1.075 3.456 1.110 -1.554
## CI_2.5  -0.578 0.728 2.905 0.796 -1.925
## CI_97.5  0.082 1.423 4.007 1.423 -1.183
## 
## $GroupPars
##         MEAN_1 COV_11
## par          0      1
## CI_2.5      NA     NA
## CI_97.5     NA     NA

###################################################
# 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)

## $overall
##    N mean_total.score sd_total.score ave.r  sd.r alpha
##  900            6.932          2.347 0.197 0.035 0.709
## 
## $itemstats
##           N  mean    sd total.r total.r_if_rm alpha_if_rm
## Item_1  900 0.570 0.495   0.604         0.443       0.673
## Item_2  900 0.689 0.463   0.518         0.351       0.690
## Item_3  900 0.526 0.500   0.531         0.352       0.691
## Item_4  900 0.892 0.310   0.422         0.305       0.698
## Item_5  900 0.748 0.435   0.522         0.367       0.687
## Item_6  900 0.543 0.498   0.569         0.398       0.682
## Item_7  900 0.759 0.428   0.530         0.379       0.685
## Item_8  900 0.803 0.398   0.516         0.375       0.686
## Item_9  900 0.777 0.417   0.489         0.337       0.692
## Item_10 900 0.626 0.484   0.548         0.378       0.685
## 
## $proportions
##             0     1
## Item_1  0.430 0.570
## Item_2  0.311 0.689
## Item_3  0.474 0.526
## Item_4  0.108 0.892
## Item_5  0.252 0.748
## Item_6  0.457 0.543
## Item_7  0.241 0.759
## Item_8  0.197 0.803
## Item_9  0.223 0.777
## Item_10 0.374 0.626

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

## 
## Call:
## lm(formula = Theta ~ covdata$group)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9871 -0.6851 -0.0427  0.7170  3.8313 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      0.04434    0.05970   0.743    0.458    
## covdata$groupg2  0.93468    0.08443  11.071   <2e-16 ***
## covdata$groupg3  1.88134    0.08443  22.284   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.034 on 897 degrees of freedom
## Multiple R-squared:  0.3563,	Adjusted R-squared:  0.3549 
## F-statistic: 248.3 on 2 and 897 DF,  p-value: < 2.2e-16

# 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)

##             AIC    SABIC       HQ      BIC    logLik      X2 df p
## rmod0  9698.483 9716.375 9718.663 9751.310 -4838.242             
## rmod1b 9472.445 9493.590 9496.294 9534.876 -4723.222 230.039  2 0

coef(rmod1a, simplify=TRUE)

## $items
##         a1      d g u
## Item_1   1 -0.466 0 1
## Item_2   1  0.193 0 1
## Item_3   1 -0.699 0 1
## Item_4   1  1.795 0 1
## Item_5   1  0.561 0 1
## Item_6   1 -0.607 0 1
## Item_7   1  0.636 0 1
## Item_8   1  0.957 0 1
## Item_9   1  0.759 0 1
## Item_10  1 -0.167 0 1
## 
## $means
## F1 
##  0 
## 
## $cov
##       F1
## F1 1.005
## 
## $lr.betas
##                F1
## (Intercept) 0.000
## groupg2     0.797
## groupg3     1.707

summary(rmod1b)

## 
## Call:
## mixedmirt(data = dat, covdata = covdata, model = 1, fixed = ~0 + 
##     items + group)
## 
## --------------
## FIXED EFFECTS:
##         Estimate Std.Error z.value
## groupg2    0.821     0.118   6.972
## groupg3    1.683     0.111  15.140
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##    F1
## F1  1

# 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)

##            AIC    SABIC       HQ      BIC    logLik      X2 df p
## mod0  9706.089 9738.620 9742.780 9802.136 -4833.044             
## mod1b 9479.361 9515.145 9519.721 9585.013 -4717.680 230.728  2 0

coef(mod1a)$lr.betas

##                    F1
## (Intercept) 0.0000000
## groupg2     0.7910307
## groupg3     1.7064184

summary(mod1b)

## 
## Call:
## mixedmirt(data = dat, covdata = covdata, model = 1, fixed = ~0 + 
##     items, itemtype = "2PL", lr.fixed = ~group)
## 
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##    F1
## F1  1
## 
## --------------
## LATENT REGRESSION FIXED EFFECTS:
## 
##                F1
## (Intercept) 0.000
## groupg2     0.767
## groupg3     1.711
## 
##             Std.Error_F1 z_F1
## (Intercept)           NA   NA
## groupg2              NaN  NaN
## groupg3              NaN  NaN

# 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]

## $GroupPars
##     MEAN_1 MEAN_2   COV_11 COV_21   COV_22
## par      0      0 1.040592      0 1.103157
## 
## $lr.betas
##                    F1        F2
## (Intercept) 0.0000000 0.0000000
## groupg2     0.7188122 0.8958095
## groupg3     1.7462013 1.6944879
## contvar     0.0431650 0.0000000

mod2b <- mixedmirt(dat, covdata, model, fixed = ~ 0 + items,
        lr.fixed = list(F1 = ~ group + contvar, F2 = ~ group))
summary(mod2b)

## 
## Call:
## mixedmirt(data = dat, covdata = covdata, model = model, fixed = ~0 + 
##     items, lr.fixed = list(F1 = ~group + contvar, F2 = ~group))
## 
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##      F1  F2
## F1 0.98 0.0
## F2 0.00 1.1
## 
## --------------
## LATENT REGRESSION FIXED EFFECTS:
## 
##                F1    F2
## (Intercept) 0.000 0.000
## groupg2     0.687 0.811
## groupg3     1.726 1.629
## contvar     0.043 0.000
## 
##             Std.Error_F1 Std.Error_F2   z_F1   z_F2
## (Intercept)           NA           NA     NA     NA
## groupg2            0.116        0.117  5.916  6.914
## groupg3            0.147        0.110 11.757 14.879
## contvar            0.052           NA  0.822     NA

####################################################
## 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)

## $overall
##     N mean_total.score sd_total.score ave.r  sd.r alpha
##  2000            5.414          2.749  0.25 0.019 0.769
## 
## $itemstats
##            N  mean    sd total.r total.r_if_rm alpha_if_rm
## Item_1  2000 0.424 0.494   0.573         0.433       0.750
## Item_2  2000 0.560 0.497   0.602         0.467       0.745
## Item_3  2000 0.373 0.484   0.547         0.405       0.754
## Item_4  2000 0.781 0.414   0.524         0.402       0.754
## Item_5  2000 0.577 0.494   0.585         0.447       0.748
## Item_6  2000 0.378 0.485   0.571         0.434       0.750
## Item_7  2000 0.598 0.491   0.568         0.428       0.751
## Item_8  2000 0.652 0.476   0.568         0.432       0.750
## Item_9  2000 0.620 0.486   0.577         0.440       0.749
## Item_10 2000 0.453 0.498   0.584         0.445       0.748
## 
## $proportions
##             0     1
## Item_1  0.577 0.424
## Item_2  0.440 0.560
## Item_3  0.627 0.373
## Item_4  0.219 0.781
## Item_5  0.424 0.577
## Item_6  0.622 0.378
## Item_7  0.402 0.598
## Item_8  0.348 0.652
## Item_9  0.380 0.620
## Item_10 0.547 0.453

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

## 
## Call:
## mixedmirt(data = dat, covdata = covdata, model = 1, fixed = ~0 + 
##     items, random = ~1 | group)
## 
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##      F1
## F1 1.02
## 
## $group
##           COV_group
## COV_group     0.853

# 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)

## 
## Call:
## mixedmirt(data = dat, covdata = covdata, model = 1, fixed = ~0 + 
##     items + group_pred, random = ~1 | group)
## 
## --------------
## FIXED EFFECTS:
##            Estimate Std.Error z.value
## group_pred    1.017      0.05  20.282
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##      F1
## F1 1.04
## 
## $group
##           COV_group
## COV_group    0.0335

anova(mod1, mod2)

##           AIC    SABIC       HQ      BIC    logLik      X2 df p
## mod1 23548.73 23577.82 23573.41 23615.94 -11762.36             
## mod2 22744.83 22776.35 22771.57 22817.65 -11359.42 805.894  1 0

# 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)

## 
## Call:
## mixedmirt(data = dat, covdata = covdata, model = 1, fixed = ~0 + 
##     items, lr.random = ~1 | group)
## 
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##      F1
## F1 1.54
## 
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $group
##           COV_group
## COV_group      1.45

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

## 
## Call:
## mixedmirt(data = dat, covdata = covdata, model = 1, fixed = ~0 + 
##     items + group_pred, lr.random = ~1 | group)
## 
## --------------
## FIXED EFFECTS:
##            Estimate Std.Error z.value
## group_pred    1.036     0.044  23.766
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##      F1
## F1 1.12
## 
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $group
##           COV_group
## COV_group       0.1

anova(mod1b, mod2b)

##            AIC    SABIC       HQ      BIC    logLik      X2 df p
## mod1b 23635.08 23661.75 23657.71 23696.69 -11806.54             
## mod2b 22749.76 22778.84 22774.44 22816.97 -11362.88 887.327  1 0

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

## 
## Call:
## mixedmirt(data = dat, covdata = covdata, model = 1, fixed = ~0 + 
##     items, itemtype = "2PL", lr.random = ~1 | group)
## 
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##    F1
## F1  1
## 
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $group
##           COV_group
## COV_group         1

anova(mod1b, mod3)

##            AIC    SABIC       HQ      BIC    logLik     X2 df p
## mod1b 23635.08 23661.75 23657.71 23696.69 -11806.54            
## mod3  23571.63 23620.11 23612.76 23683.65 -11765.82 81.452  9 0

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

##         group random_intercept
## G1  1.0631957       1.51178117
## G2 -0.4952833       0.38984324
## G3 -0.3412609      -0.62124058
## G4 -2.3189019      -2.21469989
## G5  0.7817904       1.12493092
## G6 -0.7084221      -0.04493361

## End(No test)

[Package mirt version 1.40 Index]