SimDesign

Simulation: Comparing two estimators and how effectively they recover population parameters

Item factor analysis (IFA) is a common method for determine the effectiveness and structure of items in psychological tests. Two approaches currently exist: full-information IFA via method from the item response theory framework, and limited-information IFA from the structural equation modeling framework. The question is, under which circumstances should either method be used? Here we explore a simple one factor test structure to determine which approach can recover the item parameters when varying test length and sample size.

Unfortunately, the mirt package uses a slightly different parameterization of the IFA problem, implementing a logistic model rather than a normal ogive model. Historically, the response curves have been equated by applying a scaling correction of \(D = 1.702\) to help fix this issue, which as can be seen below makes the response curves very similar (as it turns out, logistic models have slightly thicker tails than ogive models).

To give the following simulation the benefit of the doubt, the data will be generated from the normal ogive model (implying that lavaan is fitting the correct model) while mirt is only fitting an approximation to this model by translating the parameters. From this we will see whether FIML with a logistic model will better approximate the normal ogive parameters compared to the limited information DWLS approach available in lavaan.

P_logit <- function(a, d, Theta) exp(a * Theta + d) / (1 + exp(a * Theta + d))
P_ogive <- function(a, d, Theta) pnorm(a * Theta + d)

Theta <- seq(-5,5,length.out=200)
a <- 0.5
d <- -.5
D <- 1.702

example <- data.frame(Theta, logit=P_logit(a*D, d*D, Theta), ogive=P_ogive(a, d, Theta))
plot(Theta, example$logit, type = 'l', ylab = 'P', las=1)
lines(Theta, example$ogive, col = 'red')

Define the conditions

Start by defining the conditions to be studied. Because we are interested in a number of slopes and intercepts it is often convenient to save long strings of numbers to a single object to be passed to fixed_objects. Naturally, these could be included in the simulation source code directly, but often will take up a large amount of space (could be hundreds of lines), therefore passing a list object to runSimulation() may be more convenient.

library(SimDesign)
# SimFunctions()

sample_sizes <- c(250, 500, 1000)
nitems <- c(10, 20)
Design <- createDesign(sample_size = sample_sizes, 
                       nitems = nitems)

# create list of additional parameters which are fixed across conditions
set.seed(1)
pars_10 <- rbind(a = round(rlnorm(10, .3, .5)/1.702, 2),
                 d = round(rnorm(10, 0, .5)/1.702, 2))
pars_20 <- rbind(a = round(rlnorm(20, .3, .5)/1.702, 2),
                 d = round(rnorm(20, 0, .5)/1.702, 2))
(pars <- list(ten=pars_10, twenty=pars_20))

## $ten
##   [,1] [,2]  [,3]  [,4] [,5]  [,6] [,7] [,8] [,9] [,10]
## a 0.58 0.87  0.52  1.76 0.94  0.53 1.01 1.15 1.06  0.68
## d 0.44 0.11 -0.18 -0.65 0.33 -0.01 0.00 0.28 0.24  0.17
## 
## $twenty
##    [,1]  [,2] [,3] [,4]  [,5]  [,6] [,7] [,8]  [,9] [,10] [,11] [,12] [,13]
## a  1.26  1.17 0.82 0.29  1.08  0.77 0.73 0.38  0.62  0.98  1.56  0.75  0.96
## d -0.05 -0.07 0.20 0.16 -0.20 -0.21 0.11 0.23 -0.03  0.26  0.12 -0.18  0.10
##   [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## a  0.77  0.40  0.64  0.65  0.77  1.37  1.16
## d -0.33  0.42  0.58 -0.11 -0.31  0.17 -0.04

# The above slopes, when standardized, give the following factor loadings and commonalities (not run): 
a10 <- pars_10['a', ]; a20 <- pars_20['a', ]
load10 <- a10 / sqrt(1 + a10^2)
load20 <- a20 / sqrt(1 + a20^2)
rbind(F=load10, h2=load10^2)
rbind(F=load20, h2=load20^2)

Define the functions

As usual, define the functions of interest. Here we make sure that lavaan and mirt are loaded by passing them to the packages argument (required when running simulations in parallel). Here we only collect information on the slope parameters, mainly because they are the most interesting to study anyway (intercepts become more important in IRT methods).

Generate <- function(condition, fixed_objects = NULL) {

    N <- condition$sample_size
    nitems <- condition$nitems
    nitems_name <- ifelse(nitems == 10, 'ten', 'twenty')

    #extract objects from fixed_objects
    a <- fixed_objects[[nitems_name]]['a', ]
    d <- fixed_objects[[nitems_name]]['d', ]

    dat <- matrix(NA, N, nitems)
    colnames(dat) <- paste0('item_', 1:nitems)
    Theta <- rnorm(N)
    for(j in 1:nitems){
        p <- P_ogive(a[j], d[j], Theta)
        for(i in 1:N)
            dat[i,j] <- sample(c(1,0), 1, prob = c(p[i], 1 - p[i]))
    }
    as.data.frame(dat) #data.frame works nicer with lavaan
}

Analyse <- function(condition, dat, fixed_objects = NULL) {
    nitems <- condition$nitems

    # (optional) could use better starting values from fixed_objects here too
    mod <- mirt(dat, 1L, verbose=FALSE)
    if(!extract.mirt(mod, 'converged')) stop('mirt did not converge')
    cfs <- coef(mod, simplify = TRUE, digits = Inf)
    FIML_as <- cfs$items[,1L] / 1.702

    lavmod <- paste0('F =~ ', paste0('NA*', colnames(dat)[1L], ' + '),
                     paste0(colnames(dat)[-1L], collapse = ' + '),
                     '\nF ~~ 1*F')
    lmod <- sem(lavmod, dat, ordered = colnames(dat))
    if(!lavInspect(lmod, 'converged')) stop('lavaan did not converge')
    cfs2 <- coef(lmod)
    DWLS_alpha <- cfs2[1L:nitems]
    const <- sqrt(1 - DWLS_alpha^2)
    DWLS_as <- DWLS_alpha / const

    ret <- c(FIML_as=unname(FIML_as), DWLS_as=unname(DWLS_as))
    ret
}

Summarise <- function(condition, results, fixed_objects = NULL) {
    nitems <- condition$nitems
    nitems_name <- ifelse(nitems == 10, 'ten', 'twenty')

    #extract objects from fixed_objects
    a <- fixed_objects[[nitems_name]]['a', ]
    pop <- c(a, a)

    obt_bias <- bias(results, pop)
    obt_RMSE <- RMSE(results, pop)
    ret <- c(bias=obt_bias, RMSE=obt_RMSE)
    ret
}

Notice that a manual error is throw whenever lavaan or mirt objects reach their maximum number of iterations. This is required because the objects themselves do not throw error messages, but the data should be redrawn anyway to ensure that the simulation only contains models which have successfully converged.

Run the simulation

Because this simulation takes considerably longer it is recommended to pass the save = TRUE to temporarily save results in case of power outages. Results can be continued by running the identical simulation code as the initial run, and the function will automatically detect whether any temp files are available and resume the simulation at the previously saved location.

res <- runSimulation(Design, replications=100, verbose=FALSE, save=TRUE, parallel=TRUE,
                     generate=Generate, analyse=Analyse, summarise=Summarise, 
                     packages=c('mirt', 'lavaan'), fixed_objects=pars)

## This is lavaan 0.6-15
## lavaan is FREE software! Please report any bugs.

res

## # A tibble: 6 × 86
##   sample_size nitems bias.FIML_as1 bias.FIML_as2 bias.FIML_as3 bias.FIML_as4
##         <dbl>  <dbl>         <dbl>         <dbl>         <dbl>         <dbl>
## 1         250     10     -0.012443     0.029311     -0.013601      0.16573  
## 2         500     10     -0.016544     0.0032587    -0.016734      0.040801 
## 3        1000     10     -0.014192    -0.014803     -0.017341      0.093070 
## 4         250     20      0.046401     0.034873      0.0065848    -0.010772 
## 5         500     20      0.041533     0.055585      0.0065930    -0.0067927
## 6        1000     20      0.025724     0.037974     -0.0089383    -0.0048665
## # ℹ 80 more variables: bias.FIML_as5 <dbl>, bias.FIML_as6 <dbl>,
## #   bias.FIML_as7 <dbl>, bias.FIML_as8 <dbl>, bias.FIML_as9 <dbl>,
## #   bias.FIML_as10 <dbl>, bias.DWLS_as1 <dbl>, bias.DWLS_as2 <dbl>,
## #   bias.DWLS_as3 <dbl>, bias.DWLS_as4 <dbl>, bias.DWLS_as5 <dbl>,
## #   bias.DWLS_as6 <dbl>, bias.DWLS_as7 <dbl>, bias.DWLS_as8 <dbl>,
## #   bias.DWLS_as9 <dbl>, bias.DWLS_as10 <dbl>, RMSE.FIML_as1 <dbl>,
## #   RMSE.FIML_as2 <dbl>, RMSE.FIML_as3 <dbl>, RMSE.FIML_as4 <dbl>, …

Analyze the results

Sometimes, reshaping and indexing your output can be very helpful. Here we break the analysis into two parts, though other strategies are certainly possible. Because analyzing simulations is a lot like analyzing empirical data no one strategy may be the best; you have to use judgment.

For this particular analysis, we can see that the res object contains NA values for slope parameters that were not applicable. For ease of manipulating the results it often will be convenient to subset the results so that these NA’s are no longer required.

res10 <- subset(res, nitems == 10)
res10 <- res10[,!as.vector(is.na(res10[1L, ]))]
res20 <- subset(res, nitems == 20)

Ten items

# bias in slopes
names10 <- colnames(res10)
bias_as_fiml <- t(res10[,grepl('bias\\.', names10) & grepl('\\_as', names10) & 
                       grepl('FIML', names10)])
colnames(bias_as_fiml) <- sample_sizes
rownames(bias_as_fiml) <- pars_10['a', ]

bias_as_dwls <- t(res10[,grepl('bias\\.', names10) & grepl('\\_as', names10) & 
                       grepl('DWLS', names10)])
colnames(bias_as_dwls) <- sample_sizes
rownames(bias_as_dwls) <- pars_10['a', ]

(out <- list(FIML=bias_as_fiml, DWLS=bias_as_dwls))

## $FIML
##           250       500     1000
## 0.58 -0.01244 -1.65e-02 -0.01419
## 0.87  0.02931  3.26e-03 -0.01480
## 0.52 -0.01360 -1.67e-02 -0.01734
## 1.76  0.16573  4.08e-02  0.09307
## 0.94  0.02765  4.31e-05  0.00473
## 0.53 -0.01660  6.86e-04 -0.01886
## 1.01  0.02157 -1.13e-02 -0.00486
## 1.15  0.04290  1.62e-02  0.01579
## 1.06  0.01421 -2.48e-03  0.01178
## 0.68 -0.00614 -3.36e-03 -0.02256
## 
## $DWLS
##           250       500      1000
## 0.58 0.000579 -0.003796 -0.002509
## 0.87 0.043924  0.016339 -0.001135
## 0.52 0.008458  0.003640  0.000854
## 1.76 0.090144 -0.011427  0.025006
## 0.94 0.036430  0.006168  0.009948
## 0.53 0.005979  0.022785  0.002215
## 1.01 0.032245 -0.001543  0.000963
## 1.15 0.038322  0.012407  0.011271
## 1.06 0.016596  0.000158  0.012558
## 0.68 0.014988  0.014123 -0.005533

sapply(out, colMeans)

##         FIML    DWLS
## 250  0.02526 0.02877
## 500  0.00106 0.00589
## 1000 0.00327 0.00536

# RMSE in slopes
RMSE_as_fiml <- t(res10[,grepl('RMSE\\.', names10) & grepl('\\_as', names10) & 
                       grepl('FIML', names10)])
colnames(RMSE_as_fiml) <- sample_sizes
rownames(RMSE_as_fiml) <- pars_10['a', ]

RMSE_as_dwls <- t(res10[,grepl('RMSE\\.', names10) & grepl('\\_as', names10) & 
                       grepl('DWLS', names10)])
colnames(RMSE_as_dwls) <- sample_sizes
rownames(RMSE_as_dwls) <- pars_10['a', ]

(out <- list(FIML=RMSE_as_fiml, DWLS=RMSE_as_dwls))

## $FIML
##        250    500   1000
## 0.58 0.119 0.0795 0.0623
## 0.87 0.140 0.1096 0.0811
## 0.52 0.136 0.0833 0.0560
## 1.76 0.487 0.2712 0.2320
## 0.94 0.174 0.1084 0.0865
## 0.53 0.110 0.0835 0.0551
## 1.01 0.196 0.1246 0.0868
## 1.15 0.228 0.1402 0.1032
## 1.06 0.186 0.1348 0.0941
## 0.68 0.137 0.0999 0.0644
## 
## $DWLS
##        250    500   1000
## 0.58 0.117 0.0766 0.0602
## 0.87 0.141 0.1071 0.0768
## 0.52 0.138 0.0835 0.0537
## 1.76 0.424 0.2684 0.1987
## 0.94 0.171 0.1071 0.0852
## 0.53 0.108 0.0881 0.0515
## 1.01 0.190 0.1193 0.0834
## 1.15 0.216 0.1382 0.0967
## 1.06 0.177 0.1281 0.0885
## 0.68 0.137 0.1005 0.0591

sapply(out, colMeans)

##        FIML   DWLS
## 250  0.1910 0.1818
## 500  0.1235 0.1217
## 1000 0.0922 0.0854

The methods appeared to recover the slope parameters fairly well, became progressively less biased as the sample size increased, and parameters were generally recovered with greater efficiency as \(N\) increased as well. Additionally, there appeared to be an effect relating to the size of the parameters. For IRT parameters, very large slopes appeared to be recovered with greater bias compared to the DWLS, though both estimators had more difficulty recovering the more extreme slopes (see the RMSE plot below).

library(ggplot2)
plt <- data.frame(pars = c(pars_10['a', ], pars_10['a', ]), 
                  RMSE = c(RMSE_as_fiml[,'1000'], RMSE_as_dwls[,'1000']),
                  bias = c(bias_as_fiml[,'1000'], bias_as_dwls[,'1000']),
                  estimator = rep(c('FIML', 'DWLS'), each = 10))
ggplot(plt, aes(pars, bias, colour=estimator)) + geom_point(size=2) + facet_wrap(~estimator) + 
    ggtitle('slope sizes by bias for FIML and DWLS estimators')

ggplot(plt, aes(pars, RMSE, colour=estimator)) + geom_point(size=2) + facet_wrap(~estimator) + 
    ggtitle('slope sizes by RMSE for FIML and DWLS estimators')

Twenty items

# bias in slopes
names20 <- colnames(res20)
bias_as_fiml <- t(res20[,grepl('bias\\.', names20) & grepl('\\_as', names20) & 
                       grepl('FIML', names20)])
colnames(bias_as_fiml) <- sample_sizes
rownames(bias_as_fiml) <- pars_20['a', ]

bias_as_dwls <- t(res20[,grepl('bias\\.', names20) & grepl('\\_as', names20) & 
                       grepl('DWLS', names20)])
colnames(bias_as_dwls) <- sample_sizes
rownames(bias_as_dwls) <- pars_20['a', ]

(out <- list(FIML=bias_as_fiml, DWLS=bias_as_dwls))

## $FIML
##            250      500     1000
## 1.26  0.046401  0.04153  0.02572
## 1.17  0.034873  0.05558  0.03797
## 0.82  0.006585  0.00659 -0.00894
## 0.29 -0.010772 -0.00679 -0.00487
## 1.08  0.010626 -0.00121  0.01520
## 0.77  0.007025  0.00204 -0.01569
## 0.73 -0.007916 -0.01081  0.00328
## 0.38 -0.008330 -0.01382 -0.01155
## 0.62  0.002281 -0.01089 -0.00461
## 0.98  0.023817  0.01743  0.01070
## 1.56  0.117019  0.09121  0.04833
## 0.75  0.032826 -0.01756 -0.01007
## 0.96  0.000836 -0.00949  0.01265
## 0.77  0.007326  0.00661 -0.00582
## 0.4   0.006086  0.00413 -0.00196
## 0.64  0.019987 -0.00735  0.00397
## 0.65 -0.027643 -0.02360 -0.01483
## 0.77 -0.015244 -0.00893 -0.00165
## 1.37  0.056203  0.05152  0.03943
## 1.16  0.059910  0.02371  0.02496
## 
## $DWLS
##            250       500      1000
## 1.26  0.032647  0.025929  0.004138
## 1.17  0.031420  0.044931  0.022472
## 0.82  0.022136  0.019485 -0.001492
## 0.29  0.005751  0.008820  0.008835
## 1.08  0.012066 -0.002186  0.006281
## 0.77  0.023228  0.014837 -0.005545
## 0.73  0.011102  0.004688  0.015692
## 0.38  0.012220  0.002113  0.004824
## 0.62  0.020449  0.008463  0.011361
## 0.98  0.030778  0.018691  0.006861
## 1.56  0.074912  0.048274 -0.001777
## 0.75  0.049689 -0.003145  0.001264
## 0.96  0.009663 -0.003143  0.013558
## 0.77  0.016502  0.015276  0.001070
## 0.4   0.020243  0.015208  0.008523
## 0.64  0.028757  0.000215  0.006254
## 0.65 -0.006756 -0.006628  0.000775
## 0.77  0.000884  0.001665  0.007025
## 1.37  0.038533  0.026891  0.004607
## 1.16  0.052790  0.016639  0.011212

sapply(out, colMeans)

##         FIML   DWLS
## 250  0.01809 0.0244
## 500  0.00950 0.0129
## 1000 0.00711 0.0063

# RMSE in slopes
RMSE_as_fiml <- t(res20[,grepl('RMSE\\.', names20) & grepl('\\_as', names20) & 
                       grepl('FIML', names20)])
colnames(RMSE_as_fiml) <- sample_sizes
rownames(RMSE_as_fiml) <- pars_20['a', ]

RMSE_as_dwls <- t(res20[,grepl('RMSE\\.', names20) & grepl('\\_as', names20) & 
                       grepl('DWLS', names20)])
colnames(RMSE_as_dwls) <- sample_sizes
rownames(RMSE_as_dwls) <- pars_20['a', ]

(out <- list(FIML=RMSE_as_fiml, DWLS=RMSE_as_dwls))

## $FIML
##         250    500   1000
## 1.26 0.2133 0.1550 0.0989
## 1.17 0.1769 0.1466 0.1062
## 0.82 0.1378 0.0956 0.0621
## 0.29 0.0857 0.0666 0.0435
## 1.08 0.1605 0.1194 0.0846
## 0.77 0.1328 0.0896 0.0582
## 0.73 0.1207 0.0939 0.0657
## 0.38 0.1090 0.0742 0.0476
## 0.62 0.1185 0.0824 0.0561
## 0.98 0.1658 0.1134 0.0797
## 1.56 0.2527 0.2204 0.1341
## 0.75 0.1443 0.0932 0.0685
## 0.96 0.1483 0.1079 0.0831
## 0.77 0.1408 0.0862 0.0654
## 0.4  0.1239 0.0831 0.0466
## 0.64 0.1232 0.0910 0.0688
## 0.65 0.1163 0.0814 0.0648
## 0.77 0.1471 0.0867 0.0696
## 1.37 0.2284 0.1826 0.1314
## 1.16 0.2077 0.1294 0.0844
## 
## $DWLS
##         250    500   1000
## 1.26 0.1944 0.1423 0.0893
## 1.17 0.1680 0.1365 0.0961
## 0.82 0.1383 0.0930 0.0575
## 0.29 0.0893 0.0692 0.0464
## 1.08 0.1485 0.1140 0.0802
## 0.77 0.1333 0.0910 0.0557
## 0.73 0.1225 0.0904 0.0665
## 0.38 0.1131 0.0736 0.0477
## 0.62 0.1195 0.0823 0.0559
## 0.98 0.1598 0.1101 0.0746
## 1.56 0.2230 0.1879 0.1181
## 0.75 0.1456 0.0922 0.0658
## 0.96 0.1420 0.1046 0.0804
## 0.77 0.1370 0.0842 0.0638
## 0.4  0.1273 0.0859 0.0481
## 0.64 0.1265 0.0892 0.0686
## 0.65 0.1129 0.0758 0.0619
## 0.77 0.1447 0.0833 0.0690
## 1.37 0.2134 0.1678 0.1160
## 1.16 0.1943 0.1206 0.0759

sapply(out, colMeans)

##       FIML   DWLS
## 250  0.153 0.1477
## 500  0.110 0.1047
## 1000 0.076 0.0719

Again, the estimators appeared to recover the parameters with similar precision and bias. However, the effect of parameter size now is more evident in the FIML estimator. Larger slopes indeed cause more progressively bias and larger RMSE values compared to the DWLS estimator. In general, larger slopes when using FIML estimation will be under-estimated, and therefore have larger RMSEs than the DWLS approach.

library(ggplot2)
plt <- data.frame(pars = c(pars_20['a', ], pars_20['a', ]), 
                  RMSE = c(RMSE_as_fiml[,'1000'], RMSE_as_dwls[,'1000']),
                  bias = c(bias_as_fiml[,'1000'], bias_as_dwls[,'1000']),
                  estimator = rep(c('FIML', 'DWLS'), each = 20))
ggplot(plt, aes(pars, bias, colour=estimator)) + geom_point(size=2) + facet_wrap(~estimator) + 
    ggtitle('slope sizes by bias for FIML and DWLS estimators')

ggplot(plt, aes(pars, RMSE, colour=estimator)) + geom_point(size=2) + facet_wrap(~estimator) + 
    ggtitle('slope sizes by RMSE for FIML and DWLS estimators')