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 )

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

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.

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

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

# NOT RUN {
#############################################################################
# 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