# NOT RUN {
#############################################################################
# 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 <- sirt::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 <- sirt::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 <- sirt::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 <- sirt::mcmc.3pno.testlet( dat, guess.prior=1000*c(25,75), testlets=testlets,
burnin=burnin, iter=iter )
summary(mod4)
plot(mod4, ask=TRUE, layout=2)
#############################################################################
# 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( stats::rnorm( ITT, mean=0, sd=1 ), 2 )
sd0 <- 1 # sd trait
sdt <- seq( 0, 2, len=TT ) # sd testlets
# simulate theta
theta <- stats::rnorm( N, sd=sd0 )
# simulate testlets
ut <- matrix(0,nrow=N, ncol=TT )
for (tt in 1:TT){
ut[,tt] <- stats::rnorm( N, sd=sdt[tt] )
}
ut <- ut[, rep(1:TT,each=I) ]
# calculate response probability
prob <- matrix( stats::pnorm( theta + ut + matrix( b, nrow=N, ncol=ITT,
byrow=TRUE ) ), N, ITT)
Y <- (matrix( stats::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 <- sirt::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 <- sirt::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" ] )
#############################################################################
# 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( stats::rnorm( ITT, mean=0, sd=1 ), 2 )
a <- round( stats::runif( ITT, 0.5, 2 ),2)
sdt <- seq( 0, 2, len=TT ) # sd testlets
sd0 <- 1
# simulate theta
theta <- stats::rnorm( N, sd=sd0 )
# simulate testlets
ut <- matrix(0,nrow=N, ncol=TT )
for (tt in 1:TT){
ut[,tt] <- stats::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( stats::pnorm( aM*theta + ut + bM ), N, ITT)
Y <- (matrix( stats::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 <- sirt::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