#############################################################################
# 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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)
#############################################################################
# SIMULATED 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 <- rnorm( G , sd = sqrt(sig2.lev2) )
theta1 <- 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){ 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 <- pnorm( theta - b )
dat <- 1 * ( m1 > matrix( runif( n*G*I ) , n*G , I ) )
#*** estimate model
mod <- 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