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')
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)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.
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-19
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>, …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)# 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.005533sapply(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.0591sapply(out, colMeans) FIML DWLS
250 0.1910 0.1818
500 0.1235 0.1217
1000 0.0922 0.0854The 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')
# 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.011212sapply(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.0759sapply(out, colMeans) FIML DWLS
250 0.153 0.1477
500 0.110 0.1047
1000 0.076 0.0719Again, 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')