############################
## SIMULATED EXAMPLE 1: Noncompensatory and compensatory IRT models
set.seed(997)
# (1) simulate data from a two-dimensional noncompensatory
# item response model
N <- 1000 # number of persons
I <- 10 # number of items
theta0 <- rnorm( N , sd= 1 )
theta1 <- theta0 + rnorm(N , sd = .7 )
theta2 <- theta0 + rnorm(N , sd = .7 )
Q <- matrix( 1 , nrow=I,ncol=2 )
Q[ 1:(I/2) , 2 ] <- 0
Q[ I,1] <- 0
b <- matrix( rnorm( I*2 ) , I , 2 )
a <- matrix( 1 , I , 2 )
# simulate data
prob <- dat <- matrix(0 , nrow=N , ncol=I )
for (ii in 1:I){
prob[,ii] <- ( plogis( theta1 - b[ii,1] ) )^Q[ii,1]
prob[,ii] <- prob[,ii] * ( plogis( theta2 - b[ii,2] ) )^Q[ii,2]
}
dat[ prob > matrix(runif( N*I),N,I) ] <- 1
colnames(dat) <- paste0("I",1:I)
#***
# Model 1: Noncompensatory 1PL model
mod1 <- smirt(dat, Qmatrix=Q , maxiter=30 )
summary(mod1)
#***
# Model 2: Noncompensatory 2PL model
mod2 <- smirt(dat,Qmatrix=Q , est.a="2PL" , maxiter=30 )
summary(mod2)
#***
# Model 3: some fixed c and d parameters different from zero or one
c.init <- rep(0,I)
c.init[ c(3,7)] <- .2
d.init <- rep(1,I)
d.init[c(4,8)] <- .95
mod3 <- smirt( dat , Qmatrix=Q , c.init=c.init , d.init=d.init )
summary(mod3)
#***
# Model 4: some estimated c and d parameters (in parameter groups)
est.c <- c.init <- rep(0,I)
c.estpars <- c(3,6,7)
c.init[ c.estpars ] <- .2
est.c[c.estpars] <- 1
est.d <- rep(0,I)
d.init <- rep(1,I)
d.estpars <- c(6,9)
d.init[ d.estpars ] <- .95
est.d[ d.estpars ] <- d.estpars # different d parameters
mod4 <- smirt(dat,Qmatrix=Q , est.c=est.c , c.init=c.init ,
est.d=est.d , d.init=d.init )
summary(mod4)
#***
# Model 5: Unidimensional 1PL model
Qmatrix <- matrix( 1 , nrow=I , ncol=1 )
mod5 <- smirt( dat , Qmatrix=Qmatrix )
summary(mod5)
#***
# Model 6: Unidimensional 2PL model
mod6 <- smirt( dat , Qmatrix=Qmatrix , est.a="2PL" )
summary(mod6)
#***
# Model 7: Compensatory model with between item dimensionality
# between-item dimensionality
# Note that the data is simulated under the noncompensatory condition
# Therefore Model 7 should have a worse model fit than Model 1
Q1 <- Q
Q1[ 6:10 , 1] <- 0
mod7 <- smirt(dat,Qmatrix=Q1 , irtmodel="comp" , maxiter=30)
summary(mod7)
#***
# Model 8: Compensatory model with within item dimensionality
# assuming zero correlation between dimensions
variance.fixed <- as.matrix( cbind( 1,2,0) )
mod8 <- smirt(dat,Qmatrix=Q , irtmodel="comp" , variance.fixed=variance.fixed ,
maxiter=30)
summary(mod8)
#***
# Model 9: Unidimensional model with quadratic item response functions
# define theta
theta.k <- seq( - 6 , 6 , len=15 )
theta.k <- as.matrix( theta.k , ncol=1 )
# define design matrix
theta.kDES <- cbind( theta.k[,1] , theta.k[,1]^2 )
# define Q matrix
Qmatrix <- matrix( 0 , I , 2 )
Qmatrix[,1] <- 1
Qmatrix[ c(3,6,7) , 2 ] <- 1
colnames(Qmatrix) <- c("F1" , "F1sq" )
# estimate model
mod9 <- smirt(dat,Qmatrix=Qmatrix , maxiter=50 , irtmodel="comp" ,
theta.k=theta.k , theta.kDES=theta.kDES , est.a="2PL" )
summary(mod9)
#***
# Model 10: Two-dimensional item response model with latent interaction
# between dimensions
theta.k <- seq( - 6 , 6 , len=15 )
theta.k <- expand.grid( theta.k , theta.k ) # expand theta to 2 dimensions
# define design matrix
theta.kDES <- cbind( theta.k , theta.k[,1]*theta.k[,2] )
# define Q matrix
Qmatrix <- matrix( 0 , I , 3 )
Qmatrix[,1] <- 1
Qmatrix[ 6:10 , c(2,3) ] <- 1
colnames(Qmatrix) <- c("F1" , "F2" , "F1iF2" )
# estimate model
mod10 <- smirt(dat,Qmatrix=Qmatrix ,irtmodel="comp" , theta.k=theta.k ,
theta.kDES= theta.kDES , est.a="2PL" )
summary(mod10)
#****
# Model 11: Quasi Monte Carlo integration
Qmatrix <- matrix( 1 , I , 1 )
mod11 <- smirt( dat , irtmodel="comp" , Qmatrix=Qmatrix ,
qmcnodes=1000 )
summary(mod11)
Run the code above in your browser using DataLab