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.
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),
...
)a matrix or data.frame that consists of
numerically ordered data, with missing data coded as NA
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
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
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
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
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')
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.
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)
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
a list indicating parameter equality constrains. See mirt for
more detail
used for parameter starting values. See mirt for more detail
logical; return the design matrices before they have (potentially) been reassigned?
logical; compute the standard errors by approximating the information matrix using the MHRM algorithm? Default is TRUE
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
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
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.
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
# \donttest{
# 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.326 8130.107 8134.469 8171.387 -4042.663 692.217 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: 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: 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.326 8130.107 8134.469 8171.387 -4042.663
#> mod1b 7975.487 8007.270 8014.651 8077.128 -3965.743 153.84 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 7975.487 8007.270 8014.651 8077.128 -3965.743
#> mod2 8090.269 8110.494 8115.192 8154.950 -4031.134 -130.782 -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.05644 0.00341 0.04796 -0.00576 -0.04141 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : NULL
#> .. ..$ : chr "F1"
#> $ group: num [1:50, 1:2] -1.62 -1.25 -1.34 -1.55 -1.5 ...
#> ..- 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 20484.33 20525.88 20542.52 20633.82 -10208.16
#> LLTM 25464.27 25467.94 25469.41 25477.47 -12729.14 -5041.948 -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 3264.979 3276.155 3287.014 3320.577 -1618.489
#> mod2 3258.379 3271.950 3285.136 3325.891 -1612.190 12.6 3 0.006
# 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.610 9493.755 9496.459 9535.041 -4723.305 229.873 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 9481.042 9516.826 9521.402 9586.694 -4718.521 229.047 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.039167 0 1.103858
#>
#> $lr.betas
#> F1 F2
#> (Intercept) 0.00000000 0.000000
#> groupg2 0.71822234 0.906013
#> groupg3 1.74609218 1.695760
#> contvar 0.04128287 0.000000
#>
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.974 0.00
#> F2 0.000 1.12
#>
#> --------------
#> LATENT REGRESSION FIXED EFFECTS:
#>
#> F1 F2
#> (Intercept) 0.000 0.000
#> groupg2 0.675 0.803
#> groupg3 1.721 1.617
#> contvar 0.042 0.000
#>
#> Std.Error_F1 Std.Error_F2 z_F1 z_F2
#> (Intercept) NA NA NA NA
#> groupg2 NaN 0.132 NaN 6.072
#> groupg3 NaN 0.100 NaN 16.213
#> contvar 0.054 NA 0.781 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.80 23577.89 23573.48 23616.01 -11762.40
#> mod2 22744.85 22776.36 22771.59 22817.66 -11359.43 805.95 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 23637.58 23664.24 23660.20 23699.19 -11807.79
#> mod2b 22749.34 22778.43 22774.02 22816.55 -11362.67 890.24 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 23637.58 23664.24 23660.20 23699.19 -11807.79
#> mod3 23571.39 23619.87 23612.52 23683.41 -11765.69 84.189 9 0
head(cbind(randef(mod3)$group, random_intercept))
#> group random_intercept
#> G1 1.0316201 1.51178117
#> G2 -0.6199130 0.38984324
#> G3 -0.4634873 -0.62124058
#> G4 -2.4148282 -2.21469989
#> G5 0.7409736 1.12493092
#> G6 -0.5675186 -0.04493361
# }