Skip to content

Multiple booklets 1PL and PC items

tmatta edited this page Oct 17, 2017 · 3 revisions

We examined item parameter recovery under the following conditions: 2 (IRT models) x 2 (IRT R packages) x 3 (sample sizes) x 4 (test lengths) x 3 (test booklets)


  • 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
  • Two IRT R packages were evaluated: TAM (version 2.4-9) and mirt (version 1.25)
  • 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
  • Three test booklets were used: 5, 10, and 15 booklets

  • One hundred replications were used for each condition for the calibration

  • Summary of item parameter recovery:
    • TAM and mirt demonstrated a similar level of accuracy
    • For Rasch items:
      • b-parameter recovered well, with correlation ranging from 0.996 to 0.999, with bias ranging from -0.009 to 0.011, and with RMSE ranging from 0.056 to 0.339
    • For PC items:
      • b1-parameter recovered well, with correlation ranging from 0.869 to 0.996, with bias ranging from -0.029 to -0.001, and with RMSE ranging from 0.065 to 0.394
      • b2-parameter recovered well, with correlation ranging from 0.908 to 0.998, with bias ranging from -0.006 to 0.032, and with RMSE ranging from 0.062 to 0.39
    • In general:
      • Sample sizes of 5000 consistently produced the most accurate results
      • Four levels of test lengths performed very similarly
      • When number of booklets increased, recovery accuracy decreased

 

# 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
}
# 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 <- c(5, 10, 15)       #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(4375) # 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+13))
        
        # generate thetas
        theta <- rnorm(N, mean=0, sd=1)
        
        # assign items to blocks
        blocks <- lsasim::block_design(n_blocks = K, 
                                       item_parameters = item_pool, 
                                       item_block_matrix = NULL)
        
        #assign blocks to booklets
        books <- lsasim::booklet_design(item_block_assignment 
                                        = blocks$block_assignment, 
                                        book_design = NULL)
        
        #assign booklets to subjects 
        book_samp <- lsasim::booklet_sample(n_subj = N, 
                                            book_item_design = books, 
                                            book_prob = NULL)
        
        # generate item responses (Rasch and PCM)
        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 <- try(mirt::mirt(resp, 1, itemtype = 'Rasch', verbose = F))
        
        # fit Rasch and PC models using TAM package
        tam.mod <- NULL
        tam.mod <- try(TAM::tam.mml(resp, irtmodel = "PCM2", 
                                    control = list(maxiter = 200)))
        
        #------------------------------------------------------------------------------#
        # Item parameter extraction
        #------------------------------------------------------------------------------#
        if(!class(mirt.mod) == "try-error") {
          if(!class(tam.mod) == "try-error") {
            
            # 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)]
            
            #--------------------------------------------------------------------------#
            # 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)) 
            
            # combine results
            results <- rbind(results, itempars)
            
          }
        }
      }
    }
  }
}

 

  • Correlation, bias, and RMSE for item parameter recovery in mirt package

 

#--- Rasch items
mirt_recovery_1PL <- aggregate(cbind(corr_mirt_b, bias_mirt_b, RMSE_mirt_b) ~ N + I + K, 
                            data=results, mean, na.rm=TRUE)
names(mirt_recovery_1PL) <- c("N", "I", "K",
                              "corr_b", "bias_b", "RMSE_b")

round(mirt_recovery_1PL, 3)
##       N   I  K corr_b bias_b RMSE_b
## 1   500  40  5  0.989 -0.003  0.186
## 2  1000  40  5  0.995 -0.009  0.126
## 3  5000  40  5  0.999 -0.001  0.056
## 4   500  60  5  0.989 -0.004  0.190
## 5  1000  60  5  0.995 -0.006  0.126
## 6  5000  60  5  0.999 -0.002  0.058
## 7   500  80  5  0.990  0.002  0.182
## 8  1000  80  5  0.995 -0.006  0.129
## 9  5000  80  5  0.999 -0.003  0.057
## 10  500 100  5  0.989 -0.002  0.184
## 11 1000 100  5  0.995 -0.006  0.130
## 12 5000 100  5  0.999 -0.001  0.058
## 13  500  40 10  0.977  0.003  0.267
## 14 1000  40 10  0.989 -0.008  0.182
## 15 5000  40 10  0.998  0.000  0.081
## 16  500  60 10  0.979  0.007  0.260
## 17 1000  60 10  0.989 -0.004  0.186
## 18 5000  60 10  0.998 -0.003  0.081
## 19  500  80 10  0.978 -0.001  0.269
## 20 1000  80 10  0.990 -0.007  0.184
## 21 5000  80 10  0.998 -0.002  0.083
## 22  500 100 10  0.978  0.005  0.266
## 23 1000 100 10  0.989 -0.006  0.184
## 24 5000 100 10  0.998  0.001  0.081
## 25  500  40 15  0.967  0.001  0.326
## 26 1000  40 15  0.984  0.000  0.221
## 27 5000  40 15  0.997  0.001  0.101
## 28  500  60 15  0.966  0.011  0.331
## 29 1000  60 15  0.983 -0.009  0.229
## 30 5000  60 15  0.997 -0.003  0.099
## 31  500  80 15  0.967  0.001  0.339
## 32 1000  80 15  0.985 -0.004  0.224
## 33 5000  80 15  0.997 -0.002  0.103
## 34  500 100 15  0.967  0.004  0.333
## 35 1000 100 15  0.983 -0.007  0.226
## 36 5000 100 15  0.997 -0.002  0.101

 

#--- 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 + K, 
                            data=results, mean, na.rm=TRUE)
names(mirt_recovery_PC) <- c("N", "I", "K",
                             "corr_b1", "bias_b1", "RMSE_b1",
                             "corr_b2", "bias_b2", "RMSE_b2")
round(mirt_recovery_PC, 3)
##       N   I  K corr_b1 bias_b1 RMSE_b1 corr_b2 bias_b2 RMSE_b2
## 1   500  40  5   0.968  -0.005   0.216   0.979   0.001   0.209
## 2  1000  40  5   0.983  -0.009   0.154   0.990  -0.002   0.144
## 3  5000  40  5   0.996  -0.002   0.069   0.998   0.001   0.066
## 4   500  60  5   0.960  -0.001   0.214   0.974   0.003   0.209
## 5  1000  60  5   0.980  -0.008   0.148   0.987  -0.003   0.139
## 6  5000  60  5   0.996  -0.002   0.065   0.997  -0.002   0.064
## 7   500  80  5   0.955  -0.009   0.213   0.971   0.001   0.202
## 8  1000  80  5   0.977  -0.007   0.147   0.985  -0.006   0.142
## 9  5000  80  5   0.995  -0.002   0.065   0.997  -0.001   0.063
## 10  500 100  5   0.956  -0.004   0.210   0.969   0.000   0.203
## 11 1000 100  5   0.977  -0.010   0.148   0.984  -0.002   0.141
## 12 5000 100  5   0.995  -0.002   0.066   0.997  -0.001   0.062
## 13  500  40 10   0.933  -0.007   0.320   0.960   0.006   0.299
## 14 1000  40 10   0.967  -0.011   0.214   0.978  -0.001   0.214
## 15 5000  40 10   0.993  -0.001   0.098   0.996  -0.002   0.092
## 16  500  60 10   0.921  -0.016   0.312   0.947   0.006   0.295
## 17 1000  60 10   0.959  -0.012   0.213   0.973  -0.002   0.207
## 18 5000  60 10   0.992  -0.004   0.094   0.995  -0.003   0.091
## 19  500  80 10   0.911  -0.017   0.308   0.939   0.009   0.299
## 20 1000  80 10   0.955  -0.013   0.211   0.970   0.001   0.203
## 21 5000  80 10   0.990  -0.002   0.093   0.994   0.000   0.089
## 22  500 100 10   0.913  -0.017   0.309   0.939   0.004   0.290
## 23 1000 100 10   0.954  -0.015   0.214   0.969   0.000   0.201
## 24 5000 100 10   0.991  -0.002   0.093   0.993  -0.002   0.089
## 25  500  40 15   0.905  -0.023   0.384   0.935   0.030   0.387
## 26 1000  40 15   0.948  -0.012   0.274   0.968   0.010   0.257
## 27 5000  40 15   0.989  -0.005   0.120   0.993  -0.001   0.118
## 28  500  60 15   0.888  -0.029   0.378   0.921   0.015   0.377
## 29 1000  60 15   0.942  -0.015   0.262   0.959  -0.001   0.254
## 30 5000  60 15   0.987  -0.002   0.116   0.992  -0.002   0.112
## 31  500  80 15   0.872  -0.015   0.374   0.910   0.007   0.370
## 32 1000  80 15   0.931  -0.013   0.262   0.954  -0.004   0.252
## 33 5000  80 15   0.985  -0.004   0.116   0.991  -0.001   0.111
## 34  500 100 15   0.877  -0.013   0.382   0.908   0.010   0.366
## 35 1000 100 15   0.931  -0.021   0.268   0.954   0.000   0.249
## 36 5000 100 15   0.985  -0.003   0.117   0.990   0.001   0.111

 

  • Correlation, bias, and RMSE for item parameter recovery in TAM package

 

#--- Rasch items
tam_recovery_1PL <- aggregate(cbind(corr_tam_b, bias_tam_b, RMSE_tam_b) ~ N + I + K, 
                              data=results, mean, na.rm=TRUE)
names(tam_recovery_1PL) <- c("N", "I", "K",
                             "corr_b", "bias_b", "RMSE_b")
round(tam_recovery_1PL, 3)
##       N   I  K corr_b bias_b RMSE_b
## 1   500  40  5  0.989 -0.003  0.186
## 2  1000  40  5  0.995 -0.009  0.126
## 3  5000  40  5  0.999 -0.001  0.056
## 4   500  60  5  0.989 -0.005  0.190
## 5  1000  60  5  0.995 -0.006  0.126
## 6  5000  60  5  0.999 -0.002  0.058
## 7   500  80  5  0.990  0.002  0.182
## 8  1000  80  5  0.995 -0.007  0.129
## 9  5000  80  5  0.999 -0.004  0.057
## 10  500 100  5  0.989 -0.002  0.185
## 11 1000 100  5  0.995 -0.008  0.130
## 12 5000 100  5  0.999 -0.003  0.058
## 13  500  40 10  0.977  0.003  0.267
## 14 1000  40 10  0.989 -0.008  0.182
## 15 5000  40 10  0.998  0.000  0.081
## 16  500  60 10  0.979  0.007  0.260
## 17 1000  60 10  0.989 -0.004  0.186
## 18 5000  60 10  0.998 -0.003  0.081
## 19  500  80 10  0.978 -0.001  0.269
## 20 1000  80 10  0.990 -0.007  0.184
## 21 5000  80 10  0.998 -0.002  0.083
## 22  500 100 10  0.978  0.005  0.266
## 23 1000 100 10  0.989 -0.006  0.183
## 24 5000 100 10  0.998  0.001  0.081
## 25  500  40 15  0.967  0.001  0.326
## 26 1000  40 15  0.984  0.000  0.221
## 27 5000  40 15  0.997  0.001  0.101
## 28  500  60 15  0.966  0.011  0.330
## 29 1000  60 15  0.983 -0.009  0.229
## 30 5000  60 15  0.997 -0.002  0.099
## 31  500  80 15  0.967  0.001  0.339
## 32 1000  80 15  0.985 -0.004  0.224
## 33 5000  80 15  0.997 -0.001  0.103
## 34  500 100 15  0.967  0.005  0.333
## 35 1000 100 15  0.983 -0.007  0.226
## 36 5000 100 15  0.997 -0.001  0.101

 

#--- 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 + K, 
                             data=results, mean, na.rm=TRUE)
names(tam_recovery_PC) <- c("N", "I", "K",
                            "corr_b1", "bias_b1", "RMSE_b1",
                            "corr_b2", "bias_b2", "RMSE_b2")
round(tam_recovery_PC, 3)
##       N   I  K corr_b1 bias_b1 RMSE_b1 corr_b2 bias_b2 RMSE_b2
## 1   500  40  5   0.968  -0.005   0.216   0.979   0.001   0.209
## 2  1000  40  5   0.983  -0.009   0.154   0.990  -0.002   0.144
## 3  5000  40  5   0.996  -0.002   0.069   0.998   0.001   0.066
## 4   500  60  5   0.960  -0.002   0.214   0.974   0.002   0.208
## 5  1000  60  5   0.980  -0.008   0.148   0.987  -0.003   0.139
## 6  5000  60  5   0.996  -0.002   0.065   0.997  -0.002   0.064
## 7   500  80  5   0.955  -0.009   0.213   0.971   0.001   0.201
## 8  1000  80  5   0.977  -0.008   0.147   0.985  -0.006   0.142
## 9  5000  80  5   0.995  -0.003   0.065   0.997  -0.002   0.064
## 10  500 100  5   0.956  -0.004   0.211   0.969   0.000   0.203
## 11 1000 100  5   0.977  -0.012   0.149   0.984  -0.003   0.141
## 12 5000 100  5   0.995  -0.004   0.066   0.997  -0.003   0.062
## 13  500  40 10   0.933  -0.007   0.320   0.960   0.006   0.299
## 14 1000  40 10   0.967  -0.010   0.214   0.978  -0.001   0.214
## 15 5000  40 10   0.993  -0.001   0.098   0.996  -0.002   0.092
## 16  500  60 10   0.921  -0.016   0.312   0.947   0.006   0.295
## 17 1000  60 10   0.959  -0.012   0.213   0.973  -0.002   0.207
## 18 5000  60 10   0.992  -0.003   0.094   0.995  -0.002   0.091
## 19  500  80 10   0.911  -0.017   0.308   0.939   0.009   0.299
## 20 1000  80 10   0.955  -0.013   0.211   0.970   0.001   0.203
## 21 5000  80 10   0.990  -0.002   0.093   0.994   0.000   0.089
## 22  500 100 10   0.913  -0.016   0.309   0.939   0.004   0.290
## 23 1000 100 10   0.954  -0.015   0.214   0.969   0.000   0.201
## 24 5000 100 10   0.991  -0.002   0.093   0.993  -0.001   0.089
## 25  500  40 15   0.904  -0.027   0.394   0.935   0.032   0.390
## 26 1000  40 15   0.948  -0.012   0.273   0.968   0.009   0.257
## 27 5000  40 15   0.989  -0.005   0.120   0.993  -0.001   0.118
## 28  500  60 15   0.888  -0.028   0.378   0.921   0.015   0.377
## 29 1000  60 15   0.942  -0.015   0.262   0.959  -0.001   0.254
## 30 5000  60 15   0.987  -0.002   0.116   0.992  -0.002   0.112
## 31  500  80 15   0.869  -0.018   0.387   0.909   0.008   0.372
## 32 1000  80 15   0.931  -0.013   0.262   0.954  -0.004   0.252
## 33 5000  80 15   0.985  -0.004   0.116   0.991  -0.001   0.111
## 34  500 100 15   0.877  -0.012   0.382   0.908   0.010   0.366
## 35 1000 100 15   0.931  -0.020   0.268   0.954   0.000   0.249
## 36 5000 100 15   0.985  -0.002   0.117   0.990   0.001   0.111