#############################################################################
# EXAMPLE 1: Dataset Reading
#############################################################################
data(data.read)
dat <- data.read
I <- ncol(dat)
# set burnin and total number of iterations here (CHANGE THIS!)
burnin <- 200
iter <- 500
#***
# Model 1: 1PNO model
mod1 <- mcmc.3pno.testlet( dat , est.slope=FALSE , est.guess=FALSE ,
burnin=burnin, iter=iter )
summary(mod1)
plot(mod1,ask=TRUE) # plot MCMC chains in coda style
plot(mod1,ask=TRUE , layout=2) # plot MCMC output in different layout
#***
# Model 2: 3PNO model with Beta(5,17) prior for guessing parameters
mod2 <- mcmc.3pno.testlet( dat , guess.prior=c(5,17) ,
burnin=burnin, iter=iter )
summary(mod2)
#***
# Model 3: Rasch (1PNO) testlet model
testlets <- substring( colnames(dat) , 1 , 1 )
mod3 <- mcmc.3pno.testlet( dat , testlets=testlets , est.slope=FALSE ,
est.guess=FALSE , burnin=burnin, iter=iter )
summary(mod3)
#***
# Model 4: 3PNO testlet model with (almost) fixed guessing parameters .25
mod4 <- mcmc.3pno.testlet( dat , guess.prior=1000*c(25,75) , testlets=testlets ,
burnin=burnin, iter=iter )
summary(mod4)
plot(mod4, ask=TRUE, layout=2)
#############################################################################
# SIMULATED EXAMPLE 2: Simulated data according to the Rasch testlet model
#############################################################################
set.seed(678)
N <- 3000 # number of persons
I <- 4 # number of items per testlet
TT <- 3 # number of testlets
ITT <- I*TT
b <- round( rnorm( ITT , mean=0 , sd = 1 ) , 2 )
sd0 <- 1 # sd trait
sdt <- seq( 0 , 2 , len=TT ) # sd testlets
sdt <- sdt
# simulate theta
theta <- rnorm( N , sd = sd0 )
# simulate testlets
ut <- matrix(0,nrow=N , ncol=TT )
for (tt in 1:TT){ ut[,tt] <- rnorm( N , sd = sdt[tt] ) }
ut <- ut[ , rep(1:TT,each=I) ]
# calculate response probability
prob <- matrix( pnorm( theta + ut + matrix( b , nrow=N , ncol=ITT ,
byrow=TRUE ) ) , N, ITT)
Y <- (matrix( runif(N*ITT) , N , ITT) < prob )*1
colMeans(Y)
# define testlets
testlets <- rep(1:TT , each=I )
burnin <- 300
iter <- 1000
#***
# Model 1: 1PNO model (without testlet structure)
mod1 <- mcmc.3pno.testlet( dat=Y , est.slope=FALSE , est.guess=FALSE ,
burnin=burnin, iter=iter , testlets= testlets )
summary(mod1)
summ1 <- mod1$summary.mcmcobj
# compare item parameters
cbind( b , summ1[ grep("b" , summ1$parameter ) , "Mean" ] )
# Testlet standard deviations
cbind( sdt , summ1[ grep("sigma\.testlet" , summ1$parameter ) , "Mean" ] )
#***
# Model 2: 1PNO model (without testlet structure)
mod2 <- mcmc.3pno.testlet( dat=Y , est.slope=TRUE , est.guess=FALSE ,
burnin=burnin, iter=iter , testlets= testlets )
summary(mod2)
summ2 <- mod2$summary.mcmcobj
# compare item parameters
cbind( b , summ2[ grep("b\[" , summ2$parameter ) , "Mean" ] )
# item discriminations
cbind( sd0 , summ2[ grep("a\[" , summ2$parameter ) , "Mean" ] )
# Testlet standard deviations
cbind( sdt , summ2[ grep("sigma\.testlet" , summ2$parameter ) , "Mean" ] )
#############################################################################
# SIMULATED EXAMPLE 3: Simulated data according to the 2PNO testlet model
#############################################################################
set.seed(678)
N <- 3000 # number of persons
I <- 3 # number of items per testlet
TT <- 5 # number of testlets
ITT <- I*TT
b <- round( rnorm( ITT , mean=0 , sd = 1 ) , 2 )
a <- round( runif( ITT , 0.5 , 2 ) ,2)
sdt <- seq( 0 , 2 , len=TT ) # sd testlets
sdt <- sdt
# simulate theta
theta <- rnorm( N , sd = sd0 )
# simulate testlets
ut <- matrix(0,nrow=N , ncol=TT )
for (tt in 1:TT){ ut[,tt] <- rnorm( N , sd = sdt[tt] ) }
ut <- ut[ , rep(1:TT,each=I) ]
# calculate response probability
bM <- matrix( b , nrow=N , ncol=ITT , byrow=TRUE )
aM <- matrix( a , nrow=N , ncol=ITT , byrow=TRUE )
prob <- matrix( pnorm( aM*theta + ut + bM ) , N, ITT)
Y <- (matrix( runif(N*ITT) , N , ITT) < prob )*1
colMeans(Y)
# define testlets
testlets <- rep(1:TT , each=I )
burnin <- 500
iter <- 1500
#***
# Model 1: 2PNO model
mod1 <- mcmc.3pno.testlet( dat=Y , est.slope=TRUE , est.guess=FALSE ,
burnin=burnin, iter=iter , testlets= testlets )
summary(mod1)
summ1 <- mod1$summary.mcmcobj
# compare item parameters
cbind( b , summ1[ grep("b" , summ1$parameter ) , "Mean" ] )
# item discriminations
cbind( a , summ1[ grep("a\[" , summ1$parameter ) , "Mean" ] )
# Testlet standard deviations
cbind( sdt , summ1[ grep("sigma\.testlet" , summ1$parameter ) , "Mean" ] )
Run the code above in your browser using DataLab