One booklet 1PL and PC items - tmatta/lsasim GitHub Wiki
We examined item parameter recovery under the following conditions: 2 (IRT models) x 3 (IRT R packages) x 3 (sample sizes) x 4 (test lengths) x 1 (test booklet)
- Two types of IRT models were used: Rasch items and partial credit (PC) items
- Item parameters were randomly generated
- The bounds of the item difficulty parameter, b, are constrained to
b_bounds = (-2, 2)where -2 is the lowest generating value and 2 is the highest generating value
- Three IRT R packages were evaluated:
TAM(version 2.4-9),mirt(version 1.25), andltm(version 1.0-0) - Three sample sizes were used: 500, 1000, and 5000.
- Simulated samples were based on one ability level from distribution N(0, 1)
- Four test lengths were used: 40, 60, 80, and 100
- A single booklet was used.
- One hundred replications were used for each condition for the calibration
- Summary of item parameter recovery:
TAM,mirt, andltmdemonstrated a similar level of accuracy- For Rasch items:
- b-parameter recovered well, with correlation ranging from 0.996 to 1, with bias ranging from -0.016 to 0.003, and with RMSE ranging from 0.036 to 0.161
- For PC items:
- b1-parameter recovered well, with correlation ranging from 0.984 to 0.999, with bias ranging from -0.013 to -0.002, and with RMSE ranging from 0.041 to 0.160
- b2-parameter recovered well, with correlation ranging from 0.989 to 0.999, with bias ranging from -0.008 to 0.006, and with RMSE ranging from 0.040 to 0.157
- In general:
- Sample sizes of 5000 consistently produced the most accurate results
- Four levels of test lengths performed very similarly
# Load libraries
if(!require(lsasim)){
install.packages("lsasim")
library(lsasim) #version 1.0.1
}
if(!require(mirt)){
install.packages("mirt")
library(mirt) #version 1.25
}
if(!require(TAM)){
install.packages("TAM")
library(TAM) #version 2.4-9
}
if(!require(ltm)){
install.packages("ltm")
library(ltm) #version 1.0-0
}
if(!require(plyr)){
install.packages("plyr")
library(plyr) #version 1.8.4
}
# Set up conditions
N.cond <- c(500, 1000, 5000) #number of sample sizes
I.cond <- c(40, 60, 80, 100) #number of items
K.cond <- 1 #number of booklets
# Set up number of replications
reps <- 100
# Create space for outputs
results <- NULL
#==============================================================================#
# START SIMULATION
#==============================================================================#
for (N in N.cond) { #sample size
for (I in I.cond) { #number of total items
# generate item parameters for Rasch and PC models
set.seed(4373) # fix item parameters across replications
item_pool <- gen1PL_PC <- lsasim::item_gen(n_1pl = c(I/2, I/2),
thresholds = c(1, 2),
b_bounds = c(-2, 2))
for (K in K.cond) { #number of booklets
for (r in 1:reps) { #number of replications
set.seed(8088*(r+1))
# generate thetas
theta <- rnorm(N, mean=0, sd=1)
# assign items to block
block_bk1 <- lsasim::block_design(n_blocks = K,
item_parameters = item_pool)
#assign block to booklet
book_bk1 <- lsasim::booklet_design(item_block_assignment =
block_bk1$block_assignment,
book_design = matrix(K))
#assign booklet to subjects
book_samp <- lsasim::booklet_sample(n_subj = N,
book_item_design = book_bk1,
book_prob = NULL)
# generate item responses
cog <- lsasim::response_gen(subject = book_samp$subject,
item = book_samp$item,
theta = theta,
b_par = item_pool$b,
d_par = list(item_pool$d1,
item_pool$d2))
# extract item responses (excluding "subject" column)
resp <- cog[, c(1:I)]
#------------------------------------------------------------------------------#
# Item calibration
#------------------------------------------------------------------------------#
# fit Rasch and PC models using mirt package
mirt.mod <- NULL
mirt.mod <- mirt::mirt(resp, 1, itemtype = 'Rasch', verbose = F)
# fit Rasch and PC models using TAM package
tam.mod <- NULL
tam.mod <- TAM::tam.mml(resp, irtmodel = "PCM2", control = list(maxiter = 200))
# fit Rasch and PC models using ltm package
ltm.mod <- NULL
ltm.mod <- ltm::gpcm(resp, constraint = "rasch", IRT.param=T,
control = list(iter.qN = 1000))
#------------------------------------------------------------------------------#
# Item parameter extraction
#------------------------------------------------------------------------------#
# extract b, b1, b2 in mirt package
#--- Rasch items
mirt_b <- coef(mirt.mod, IRTpars = TRUE, simplify=TRUE)$items[ c(1:(I/2)),"b"]
#--- PC items
mirt_b1 <- coef(mirt.mod, IRTpars = TRUE, simplify=TRUE)$items[ c( (I/2+1) : I),"b1"]
mirt_b2 <- coef(mirt.mod, IRTpars = TRUE, simplify=TRUE)$items[ c( (I/2+1) : I),"b2"]
# convert TAM output into PCM parametrization
#--- Rasch items
tam_b <- tam.mod$item$xsi.item[ c(1: (I/2))]
#--- PC items
tam_b1 <- tam.mod$item$AXsi_.Cat1[ c( (I/2+1) : I)]
tam_b2 <- ((tam.mod$item$AXsi_.Cat2) - (tam.mod$item$AXsi_.Cat1))[ c( (I/2+1) : I)]
# extract Catgr.1 and Catgr.2 in ltm package
ltm.mod.par <- plyr::ldply(coef(ltm.mod), rbind)
#--- Rasch items
ltm_b <- ltm.mod.par$Catgr.1[ c(1: (I/2))]
#--- PC items
ltm_b1 <- ltm.mod.par$Catgr.1[ c( (I/2+1) : I)]
ltm_b2 <- ltm.mod.par$Catgr.2[ c( (I/2+1) : I)]
#------------------------------------------------------------------------------#
# Item parameter recovery
#------------------------------------------------------------------------------#
# summarize results
itempars <- data.frame(matrix(c(N, I, K, r), nrow=1))
colnames(itempars) <- c("N", "I", "K", "rep")
# retrieve generated item parameters
genPC.b1 <- item_pool$b + item_pool$d1
genPC.b2 <- item_pool$b + item_pool$d2
# calculate corr, bias, RMSE for item parameters in mirt pacakge
itempars$corr_mirt_b <- cor( item_pool$b[ c(1: (I/2))], mirt_b)
itempars$bias_mirt_b <- mean( mirt_b - item_pool$b[ c(1: (I/2))] )
itempars$RMSE_mirt_b <- sqrt(mean( ( mirt_b - item_pool$b[ c(1: (I/2))])^2 ))
itempars$corr_mirt_b1 <- cor( item_pool$b1[ c( (I/2+1) : I)], mirt_b1)
itempars$bias_mirt_b1 <- mean( mirt_b1 - item_pool$b1[ c( (I/2+1) : I)] )
itempars$RMSE_mirt_b1 <- sqrt(mean( ( mirt_b1 - item_pool$b1[ c( (I/2+1) : I)])^2 ))
itempars$corr_mirt_b2 <- cor( item_pool$b2[ c( (I/2+1) : I)], mirt_b2)
itempars$bias_mirt_b2 <- mean( mirt_b2 - item_pool$b2[ c( (I/2+1) : I)] )
itempars$RMSE_mirt_b2 <- sqrt(mean( ( mirt_b2 - item_pool$b2[ c( (I/2+1) : I)])^2 ))
# calculate corr, bias, RMSE for item parameters in TAM pacakge
itempars$corr_tam_b <- cor( item_pool$b[ c(1: (I/2))], tam_b)
itempars$bias_tam_b <- mean( tam_b - item_pool$b[ c(1: (I/2))] )
itempars$RMSE_tam_b <- sqrt(mean( ( tam_b - item_pool$b[ c(1:(I/2))])^2 ))
itempars$corr_tam_b1 <- cor( item_pool$b1[ c( (I/2+1) : I)], tam_b1)
itempars$bias_tam_b1 <- mean( tam_b1 - item_pool$b1[ c( (I/2+1) : I)] )
itempars$RMSE_tam_b1 <- sqrt(mean( ( tam_b1 - item_pool$b1[ c( (I/2+1) : I)])^2 ))
itempars$corr_tam_b2 <- cor( item_pool$b2[ c( (I/2+1) : I)], tam_b2)
itempars$bias_tam_b2 <- mean( tam_b2 - item_pool$b2[ c( (I/2+1) : I)] )
itempars$RMSE_tam_b2 <- sqrt(mean( ( tam_b2 - item_pool$b2[ c( (I/2+1) : I)])^2 ))
# calculate corr, bias, RMSE for item parameters in ltm pacakge
itempars$corr_ltm_b <- cor( item_pool$b[ c(1: (I/2))], ltm_b)
itempars$bias_ltm_b <- mean( ltm_b - item_pool$b[ c(1:(I/2))] )
itempars$RMSE_ltm_b <- sqrt(mean( ( ltm_b - item_pool$b[ c(1:(I/2))])^2 ))
itempars$corr_ltm_b1 <- cor( item_pool$b1[ c( (I/2+1) : I)], ltm_b1)
itempars$bias_ltm_b1 <- mean( ltm_b1 - item_pool$b1[ c( (I/2+1) : I)] )
itempars$RMSE_ltm_b1 <- sqrt(mean( ( ltm_b1 - item_pool$b1[ c( (I/2+1) : I)])^2 ))
itempars$corr_ltm_b2 <- cor( item_pool$b2[ c( (I/2+1) : I)], ltm_b2)
itempars$bias_ltm_b2 <- mean( ltm_b2 - item_pool$b2[ c( (I/2+1) : I)] )
itempars$RMSE_ltm_b2 <- sqrt(mean( ( ltm_b2 - item_pool$b2[ c( (I/2+1) : I)])^2 ))
# combine results
results <- rbind(results, itempars)
}
}
}
}
- Correlation, bias, and RMSE for item parameter recovery in
mirtpackage
#--- Rasch items
mirt_recovery_1PL <- aggregate(cbind(corr_mirt_b, bias_mirt_b, RMSE_mirt_b) ~ N + I,
data=results, mean, na.rm=TRUE)
names(mirt_recovery_1PL) <- c("Sample Size", "Test Length",
"corr_b", "bias_b", "RMSE_b")
round(mirt_recovery_1PL, 3)
## Sample Size Test Length corr_b bias_b RMSE_b
## 1 500 40 0.996 0.001 0.115
## 2 1000 40 0.998 -0.008 0.078
## 3 5000 40 1.000 -0.002 0.036
## 4 500 60 0.996 -0.004 0.118
## 5 1000 60 0.998 -0.006 0.082
## 6 5000 60 1.000 -0.001 0.037
## 7 500 80 0.996 -0.006 0.117
## 8 1000 80 0.998 -0.006 0.084
## 9 5000 80 1.000 -0.003 0.038
## 10 500 100 0.996 -0.004 0.120
## 11 1000 100 0.998 -0.008 0.083
## 12 5000 100 1.000 -0.002 0.037
#--- PC items
mirt_recovery_PC <- aggregate(cbind(corr_mirt_b1, bias_mirt_b1, RMSE_mirt_b1,
corr_mirt_b2, bias_mirt_b2, RMSE_mirt_b2) ~ N + I,
data=results, mean, na.rm=TRUE)
names(mirt_recovery_PC) <- c("Sample Size", "Test Length",
"corr_b1", "bias_b1", "RMSE_b1",
"corr_b2", "bias_b2", "RMSE_b2")
round(mirt_recovery_PC, 3)
## Sample Size Test Length corr_b1 bias_b1 RMSE_b1 corr_b2 bias_b2 RMSE_b2
## 1 500 40 0.990 -0.008 0.131 0.992 0.002 0.136
## 2 1000 40 0.995 -0.003 0.091 0.996 -0.004 0.100
## 3 5000 40 0.999 -0.002 0.042 0.999 -0.001 0.043
## 4 500 60 0.987 -0.008 0.138 0.990 -0.002 0.140
## 5 1000 60 0.994 -0.009 0.095 0.995 -0.007 0.094
## 6 5000 60 0.999 -0.002 0.042 0.999 -0.002 0.043
## 7 500 80 0.986 -0.010 0.131 0.989 0.001 0.135
## 8 1000 80 0.993 -0.006 0.092 0.994 -0.006 0.094
## 9 5000 80 0.999 -0.002 0.041 0.999 -0.001 0.042
## 10 500 100 0.984 -0.008 0.135 0.989 -0.004 0.131
## 11 1000 100 0.992 -0.009 0.095 0.994 -0.005 0.093
## 12 5000 100 0.998 -0.004 0.042 0.999 -0.003 0.040
- Correlation, bias, and RMSE for item parameter recovery in
TAMpackage
#--- Rasch items
tam_recovery_1PL <- aggregate(cbind(corr_tam_b, bias_tam_b, RMSE_tam_b) ~ N + I,
data=results, mean, na.rm=TRUE)
names(tam_recovery_1PL) <- c("Sample Size", "Test Length",
"corr_b", "bias_b", "RMSE_b")
round(tam_recovery_1PL, 3)
## Sample Size Test Length corr_b bias_b RMSE_b
## 1 500 40 0.996 0.001 0.116
## 2 1000 40 0.998 -0.007 0.078
## 3 5000 40 1.000 -0.002 0.036
## 4 500 60 0.996 0.000 0.124
## 5 1000 60 0.998 -0.004 0.085
## 6 5000 60 1.000 -0.002 0.038
## 7 500 80 0.996 -0.004 0.135
## 8 1000 80 0.998 -0.010 0.095
## 9 5000 80 1.000 -0.004 0.042
## 10 500 100 0.996 -0.006 0.148
## 11 1000 100 0.998 0.001 0.111
## 12 5000 100 1.000 0.000 0.048
#--- PC items
tam_recovery_PC <- aggregate(cbind(corr_tam_b1, bias_tam_b1, RMSE_tam_b1,
corr_tam_b2, bias_tam_b2, RMSE_tam_b2) ~ N + I,
data=results, mean, na.rm=TRUE)
names(tam_recovery_PC) <- c("Sample Size", "Test Length",
"corr_b1", "bias_b1", "RMSE_b1",
"corr_b2", "bias_b2", "RMSE_b2")
round(tam_recovery_PC, 3)
## Sample Size Test Length corr_b1 bias_b1 RMSE_b1 corr_b2 bias_b2 RMSE_b2
## 1 500 40 0.990 -0.008 0.132 0.992 0.002 0.137
## 2 1000 40 0.995 -0.003 0.091 0.996 -0.003 0.100
## 3 5000 40 0.999 -0.002 0.042 0.999 -0.001 0.043
## 4 500 60 0.987 -0.006 0.145 0.990 0.003 0.146
## 5 1000 60 0.994 -0.007 0.097 0.995 -0.004 0.097
## 6 5000 60 0.999 -0.003 0.043 0.999 -0.001 0.044
## 7 500 80 0.986 -0.010 0.145 0.989 0.005 0.150
## 8 1000 80 0.993 -0.011 0.102 0.994 -0.008 0.103
## 9 5000 80 0.999 -0.005 0.045 0.999 0.000 0.046
## 10 500 100 0.984 -0.013 0.160 0.989 -0.003 0.157
## 11 1000 100 0.992 -0.003 0.122 0.994 0.006 0.119
## 12 5000 100 0.998 -0.005 0.052 0.999 0.003 0.051
- Correlation, bias, and RMSE for item parameter recovery in
ltmpackage
#--- Rasch items
ltm_recovery_1PL <- aggregate(cbind(corr_ltm_b, bias_ltm_b, RMSE_ltm_b) ~ N + I,
data=results, mean, na.rm=TRUE)
names(ltm_recovery_1PL) <- c("Sample Size", "Test Length",
"corr_b", "bias_b", "RMSE_b")
round(ltm_recovery_1PL, 3)
## Sample Size Test Length corr_b bias_b RMSE_b
## 1 500 40 0.996 -0.002 0.118
## 2 1000 40 0.998 -0.007 0.079
## 3 5000 40 1.000 -0.004 0.037
## 4 500 60 0.996 -0.003 0.132
## 5 1000 60 0.998 -0.003 0.092
## 6 5000 60 1.000 -0.001 0.040
## 7 500 80 0.996 -0.013 0.144
## 8 1000 80 0.998 -0.007 0.109
## 9 5000 80 1.000 -0.002 0.049
## 10 500 100 0.996 -0.016 0.161
## 11 1000 100 0.998 0.003 0.117
## 12 5000 100 1.000 0.000 0.056
#--- PC items
ltm_recovery_PC <- aggregate(cbind(corr_ltm_b1, bias_ltm_b1, RMSE_ltm_b1,
corr_ltm_b2, bias_ltm_b2, RMSE_ltm_b2) ~ N + I,
data=results, mean, na.rm=TRUE)
names(ltm_recovery_PC) <- c("Sample Size", "Test Length",
"corr_b1", "bias_b1", "RMSE_b1",
"corr_b2", "bias_b2", "RMSE_b2")
round(ltm_recovery_PC, 3)
## Sample Size Test Length corr_b1 bias_b1 RMSE_b1 corr_b2 bias_b2 RMSE_b2
## 1 500 40 0.990 -0.010 0.133 0.992 -0.001 0.139
## 2 1000 40 0.995 -0.003 0.091 0.996 -0.003 0.100
## 3 5000 40 0.999 -0.004 0.042 0.999 -0.003 0.044
## 4 500 60 0.987 -0.009 0.152 0.990 0.000 0.153
## 5 1000 60 0.994 -0.008 0.103 0.995 -0.002 0.102
## 6 5000 60 0.999 -0.004 0.045 0.999 0.001 0.046
## 7 500 80 0.986 -0.020 0.155 0.989 -0.002 0.159
## 8 1000 80 0.993 -0.010 0.116 0.994 -0.003 0.115
## 9 5000 80 0.999 -0.005 0.052 0.999 0.004 0.053
## 10 500 100 0.984 -0.026 0.172 0.989 -0.012 0.169
## 11 1000 100 0.992 -0.004 0.128 0.994 0.010 0.125
## 12 5000 100 0.998 -0.008 0.060 0.999 0.005 0.059