mcmc.2pno.ml {sirt}

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

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

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

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.

See Also

S3 methods: summary.mcmc.sirt, plot.mcmc.sirt

For MCMC estimation of three-parameter (testlet) models see mcmc.3pno.testlet.

See also the MLIRT package (http://www.jean-paulfox.com).

For more flexible estimation of multilevel IRT models see the MCMCglmm and lme4 packages.

Examples

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

## End(Not run)