mcmc.2pno.ml: Random Item Response Model / Multilevel IRT Model

Description

This function enables the estimation of random item models and multilevel (or hierarchical) IRT models (Chaimongkol, Huffer & Kamata, 2007; Fox & Verhagen, 2010; van den Noortgate, de Boeck & Meulders, 2003; Asparouhov & Muthen, 2012; Muthen & Asparouhov, 2013, 2014). Dichotomous response data is supported using a probit link. Normally distributed responses can also be analyzed. See Details for a description of the implemented item response models.

Usage

mcmc.2pno.ml(dat, group, link="logit", est.b.M="h", est.b.Var="n",
    est.a.M="f", est.a.Var="n", burnin=500, iter=1000,
    N.sampvalues=1000, progress.iter=50, prior.sigma2=c(1, 0.4),
    prior.sigma.b=c(1, 1), prior.sigma.a=c(1, 1), prior.omega.b=c(1, 1),
    prior.omega.a=c(1, 0.4), sigma.b.init=.3 )

Value

A list of class mcmc.sirt with following entries:

mcmcobj: Object of class mcmc.list
summary.mcmcobj: Summary of the mcmcobj object. In this summary the Rhat statistic and the mode estimate MAP is included. The variable PercSEratio indicates the proportion of the Monte Carlo standard error in relation to the total standard deviation of the posterior distribution.
ic: Information criteria (DIC)
burnin: Number of burnin iterations
iter: Total number of iterations
theta.chain: Sampled values of $\theta_{pj}$ parameters
theta.chain: Sampled values of $u_{j}$ parameters
deviance.chain: Sampled values of Deviance values
EAP.rel: EAP reliability
person: Data frame with EAP person parameter estimates for $\theta_pj$ and their corresponding posterior standard deviations
dat: Used data frame
...: Further values

Arguments

dat: Data frame with item responses.
group: Vector of group identifiers (e.g. classes, schools or countries)
link: Link function. Choices are "logit" for dichotomous data and "normal" for data under normal distribution assumptions
est.b.M: Estimation type of $b_i$ parameters:
n - non-hierarchical prior distribution, i.e. $\omega_b$ is set to a very high value and is not estimated
h - hierarchical prior distribution with estimated distribution parameters $\mu_b$ and $\omega_b$
est.b.Var: Estimation type of standard deviations of item difficulties $b_i$.
n -- no estimation of the item variance, i.e. $\sigma_{b,i}$ is assumed to be zero
i -- item-specific standard deviation of item difficulties
j -- a joint standard deviation of all item difficulties is estimated, i.e. $\sigma_{b,1}=\ldots=\sigma_{b,I}=\sigma_b$
est.a.M: Estimation type of $a_i$ parameters:
f - no estimation of item slopes, i.e all item slopes $a_i$ are fixed at one
n - non-hierarchical prior distribution, i.e. $\omega_a=0$
h - hierarchical prior distribution with estimated distribution parameter $\omega_a$
est.a.Var: Estimation type of standard deviations of item slopes $a_i$.
n -- no estimation of the item variance
i -- item-specific standard deviation of item slopes
j -- a joint standard deviation of all item slopes is estimated, i.e. $\sigma_{a,1}=\ldots=\sigma_{a,I}=\sigma_a$
burnin: Number of burnin iterations
iter: Total number of iterations
N.sampvalues: Maximum number of sampled values to save
progress.iter: Display progress every progress.iter-th iteration. If no progress display is wanted, then choose progress.iter larger than iter.
prior.sigma2: Prior for Level 2 standard deviation $\sigma_{L2}$
prior.sigma.b: Priors for item difficulty standard deviations $\sigma_{b,i}$
prior.sigma.a: Priors for item difficulty standard deviations $\sigma_{a,i}$
prior.omega.b: Prior for $\omega_b$
prior.omega.a: Prior for $\omega_a$
sigma.b.init: Initial standard deviation for $\sigma_{b,i}$ parameters

Details

For dichotomous item responses (link="logit") of persons $p$ in group $j$ on item $i$, the probability of a correct response is defined as $$P( X_{pji}=1 | \theta_{pj} )=\Phi ( a_{ij} \theta_{pj} - b_{ij} )$$ The ability $\theta_{pj}$ is decomposed into a Level 1 and a Level 2 effect $$\theta_{pj}=u_j + e_{pj} \quad, \quad u_j \sim N ( 0, \sigma_{L2}^2 ) \quad, \quad e_{pj} \sim N ( 0, \sigma_{L1}^2 ) $$ In a multilevel IRT model (or a random item model), item parameters are allowed to vary across groups: $$ b_{ij} \sim N( b_i, \sigma^2_{b,i} ) \quad, \quad a_{ij} \sim N( a_i, \sigma^2_{a,i} ) $$ In a hierarchical IRT model, a hierarchical distribution of the (main) item parameters is assumed $$ b_{i} \sim N( \mu_b, \omega^2_{b} ) \quad, \quad a_{i} \sim N( 1, \omega^2_{a} ) $$ Note that for identification purposes, the mean of all item slopes $a_i$ is set to one. Using the arguments est.b.M, est.b.Var, est.a.M and est.a.Var defines which variance components should be estimated.

For normally distributed item responses (link="normal"), the model equations remain the same except the item response model which is now written as $$ X_{pji}=a_{ij} \theta_{pj} - b_{ij} + \varepsilon_{pji} \quad, \quad \varepsilon_{pji} \sim N( 0, \sigma^2_{res,i} ) $$

References

Asparouhov, T. & Muthen, B. (2012). General random effect latent variable modeling: Random subjects, items, contexts, and parameters. http://www.statmodel.com/papers_date.shtml.

Chaimongkol, S., Huffer, F. W., & Kamata, A. (2007). An explanatory differential item functioning (DIF) model by the WinBUGS 1.4. Songklanakarin Journal of Science and Technology, 29, 449-458.

Fox, J.-P., & Verhagen, A.-J. (2010). Random item effects modeling for cross-national survey data. In E. Davidov, P. Schmidt, & J. Billiet (Eds.), Cross-cultural Analysis: Methods and Applications (pp. 467-488), London: Routledge Academic.

Muthen, B. & Asparouhov, T. (2013). New methods for the study of measurement invariance with many groups. http://www.statmodel.com/papers_date.shtml

Muthen, B. & Asparouhov, T. (2014). Item response modeling in Mplus: A multi-dimensional, multi-level, and multi-timepoint example. In W. Linden & R. Hambleton (2014). Handbook of item response theory: Models, statistical tools, and applications. http://www.statmodel.com/papers_date.shtml

van den Noortgate, W., De Boeck, P., & Meulders, M. (2003). Cross-classification multilevel logistic models in psychometrics. Journal of Educational and Behavioral Statistics, 28, 369-386.

Examples

Run this code

if (FALSE) {
#############################################################################
# EXAMPLE 1: Dataset Multilevel data.ml1 - dichotomous items
#############################################################################
data(data.ml1)
dat <- data.ml1[,-1]
group <- data.ml1$group
# just for a try use a very small number of iterations
burnin <- 50 ; iter <- 100

#***
# Model 1: 1PNO with no cluster item effects
mod1 <- sirt::mcmc.2pno.ml( dat, group, est.b.Var="n", burnin=burnin, iter=iter )
summary(mod1)    # summary
plot(mod1,layout=2,ask=TRUE) # plot results
# write results to coda file
mcmclist2coda( mod1$mcmcobj, name="data.ml1_mod1" )

#***
# Model 2: 1PNO with cluster item effects of item difficulties
mod2 <- sirt::mcmc.2pno.ml( dat, group, est.b.Var="i", burnin=burnin, iter=iter )
summary(mod2)
plot(mod2, ask=TRUE, layout=2 )

#***
# Model 3: 2PNO with cluster item effects of item difficulties but
#          joint item slopes
mod3 <- sirt::mcmc.2pno.ml( dat, group, est.b.Var="i", est.a.M="h",
              burnin=burnin, iter=iter )
summary(mod3)

#***
# Model 4: 2PNO with cluster item effects of item difficulties and
#          cluster item effects with a jointly estimated SD
mod4 <- sirt::mcmc.2pno.ml( dat, group, est.b.Var="i", est.a.M="h",
                est.a.Var="j", burnin=burnin, iter=iter )
summary(mod4)

#############################################################################
# EXAMPLE 2: Dataset Multilevel data.ml2 - polytomous items
#            assuming a normal distribution for polytomous items
#############################################################################
data(data.ml2)
dat <- data.ml2[,-1]
group <- data.ml2$group
# set iterations for all examples (too few!!)
burnin <- 100 ; iter <- 500

#***
# Model 1: no intercept variance, no slopes
mod1 <- sirt::mcmc.2pno.ml( dat=dat, group=group, est.b.Var="n",
             burnin=burnin, iter=iter, link="normal",  progress.iter=20  )
summary(mod1)

#***
# Model 2a: itemwise intercept variance, no slopes
mod2a <- sirt::mcmc.2pno.ml( dat=dat, group=group, est.b.Var="i",
            burnin=burnin, iter=iter,link="normal",  progress.iter=20  )
summary(mod2a)

#***
# Model 2b: homogeneous intercept variance, no slopes
mod2b <- sirt::mcmc.2pno.ml( dat=dat, group=group, est.b.Var="j",
              burnin=burnin, iter=iter,link="normal",  progress.iter=20  )
summary(mod2b)

#***
# Model 3: intercept variance and slope variances
#          hierarchical item and slope parameters
mod3 <- sirt::mcmc.2pno.ml( dat=dat, group=group,
               est.b.M="h", est.b.Var="i", est.a.M="h", est.a.Var="i",
               burnin=burnin, iter=iter,link="normal",  progress.iter=20  )
summary(mod3)

#############################################################################
# EXAMPLE 3: Simulated random effects model | dichotomous items
#############################################################################
set.seed(7698)

#*** model parameters
sig2.lev2 <- .3^2   # theta level 2 variance
sig2.lev1 <- .8^2   # theta level 1 variance
G <- 100            # number of groups
n <- 20             # number of persons within a group
I <- 12             # number of items
#*** simuate theta
theta2 <- stats::rnorm( G, sd=sqrt(sig2.lev2) )
theta1 <- stats::rnorm( n*G, sd=sqrt(sig2.lev1) )
theta  <- theta1 + rep( theta2, each=n )
#*** item difficulties
b <- seq( -2, 2, len=I )
#*** define group identifier
group <- 1000 + rep(1:G, each=n )
#*** SD of group specific difficulties for items 3 and 5
sigma.item <- rep(0,I)
sigma.item[c(3,5)] <- 1
#*** simulate group specific item difficulties
b.class <- sapply( sigma.item, FUN=function(sii){ stats::rnorm( G, sd=sii ) } )
b.class <- b.class[ rep( 1:G,each=n ), ]
b <- matrix( b, n*G, I, byrow=TRUE ) + b.class
#*** simulate item responses
m1 <- stats::pnorm( theta - b )
dat <- 1 * ( m1 > matrix( stats::runif( n*G*I ), n*G, I ) )

#*** estimate model
mod <- sirt::mcmc.2pno.ml( dat, group=group, burnin=burnin, iter=iter,
            est.b.M="n", est.b.Var="i", progress.iter=20)
summary(mod)
plot(mod, layout=2, ask=TRUE )
}

Run the code above in your browser using DataLab