This function computes updated parameter and standard error estimates using multiple imputation methodology. Given a set of parameter estimates and their associated standard errors the function returns the weighted average of the overall between and within variability due to the multiple imputations according to Rubin's (1987) methodology.
Value
returns a list or data.frame containing the updated averaged parameter estimates, standard errors, and t-values with the associated degrees of freedom and two tailed p-values
References
Chalmers, R., P. (2012). mirt: A Multidimensional Item Response Theory Package for the R Environment. Journal of Statistical Software, 48(6), 1-29. doi:10.18637/jss.v048.i06
Rubin, D.B. (1987) Multiple Imputation for Nonresponse in Surveys. Wiley & Sons, New York.
Author
Phil Chalmers rphilip.chalmers@gmail.com
Examples
# \donttest{
# simulate data
set.seed(1234)
N <- 1000
# covariates
X1 <- rnorm(N); X2 <- rnorm(N)
covdata <- data.frame(X1, X2)
Theta <- matrix(0.5 * X1 + -1 * X2 + rnorm(N, sd = 0.5))
# items and response data
a <- matrix(1, 20); d <- matrix(rnorm(20))
dat <- simdata(a, d, 1000, itemtype = '2PL', Theta=Theta)
mod1 <- mirt(dat, 1, 'Rasch', covdata=covdata, formula = ~ X1 + X2)
#>
Iteration: 1, Log-Lik: -11123.874, Max-Change: 0.77434
Iteration: 2, Log-Lik: -10724.065, Max-Change: 0.53447
Iteration: 3, Log-Lik: -10569.146, Max-Change: 0.20072
Iteration: 4, Log-Lik: -10525.011, Max-Change: 0.08901
Iteration: 5, Log-Lik: -10508.664, Max-Change: 0.04954
Iteration: 6, Log-Lik: -10501.108, Max-Change: 0.03132
Iteration: 7, Log-Lik: -10497.152, Max-Change: 0.02140
Iteration: 8, Log-Lik: -10494.904, Max-Change: 0.01541
Iteration: 9, Log-Lik: -10493.551, Max-Change: 0.01150
Iteration: 10, Log-Lik: -10492.672, Max-Change: 0.00901
Iteration: 11, Log-Lik: -10492.127, Max-Change: 0.00697
Iteration: 12, Log-Lik: -10491.775, Max-Change: 0.00549
Iteration: 13, Log-Lik: -10491.541, Max-Change: 0.00449
Iteration: 14, Log-Lik: -10491.376, Max-Change: 0.00359
Iteration: 15, Log-Lik: -10491.270, Max-Change: 0.00290
Iteration: 16, Log-Lik: -10491.197, Max-Change: 0.00240
Iteration: 17, Log-Lik: -10491.146, Max-Change: 0.00196
Iteration: 18, Log-Lik: -10491.112, Max-Change: 0.00161
Iteration: 19, Log-Lik: -10491.088, Max-Change: 0.00137
Iteration: 20, Log-Lik: -10491.070, Max-Change: 0.00112
Iteration: 21, Log-Lik: -10491.058, Max-Change: 0.00093
Iteration: 22, Log-Lik: -10491.050, Max-Change: 0.00078
Iteration: 23, Log-Lik: -10491.044, Max-Change: 0.00065
Iteration: 24, Log-Lik: -10491.040, Max-Change: 0.00054
Iteration: 25, Log-Lik: -10491.037, Max-Change: 0.00046
Iteration: 26, Log-Lik: -10491.035, Max-Change: 0.00039
Iteration: 27, Log-Lik: -10491.034, Max-Change: 0.00032
Iteration: 28, Log-Lik: -10491.033, Max-Change: 0.00027
Iteration: 29, Log-Lik: -10491.032, Max-Change: 0.00023
Iteration: 30, Log-Lik: -10491.032, Max-Change: 0.00019
Iteration: 31, Log-Lik: -10491.031, Max-Change: 0.00026
Iteration: 32, Log-Lik: -10491.031, Max-Change: 0.00014
Iteration: 33, Log-Lik: -10491.031, Max-Change: 0.00012
Iteration: 34, Log-Lik: -10491.031, Max-Change: 0.00010
coef(mod1, simplify=TRUE)
#> $items
#> a1 d g u
#> Item_1 1 -0.409 0 1
#> Item_2 1 0.491 0 1
#> Item_3 1 0.313 0 1
#> Item_4 1 1.965 0 1
#> Item_5 1 1.753 0 1
#> Item_6 1 -0.246 0 1
#> Item_7 1 -1.077 0 1
#> Item_8 1 0.533 0 1
#> Item_9 1 -1.232 0 1
#> Item_10 1 0.603 0 1
#> Item_11 1 -0.404 0 1
#> Item_12 1 1.238 0 1
#> Item_13 1 1.033 0 1
#> Item_14 1 1.524 0 1
#> Item_15 1 -0.548 0 1
#> Item_16 1 2.075 0 1
#> Item_17 1 -0.695 0 1
#> Item_18 1 -1.200 0 1
#> Item_19 1 0.121 0 1
#> Item_20 1 0.523 0 1
#>
#> $means
#> F1
#> 0
#>
#> $cov
#> F1
#> F1 0.215
#>
#> $lr.betas
#> F1
#> (Intercept) 0.000
#> X1 0.527
#> X2 -1.036
#>
# draw plausible values for secondary analyses
pv <- fscores(mod1, plausible.draws = 10)
pvmods <- lapply(pv, function(x, covdata) lm(x ~ covdata$X1 + covdata$X2),
covdata=covdata)
# compute Rubin's multiple imputation average
so <- lapply(pvmods, summary)
par <- lapply(so, function(x) x$coefficients[, 'Estimate'])
SEpar <- lapply(so, function(x) x$coefficients[, 'Std. Error'])
averageMI(par, SEpar)
#> par SEpar t df p
#> (Intercept) 0.003 0.016 0.209 198.552 0.209
#> covdata$X1 0.528 0.018 28.554 63.384 0
#> covdata$X2 -1.037 0.021 -49.311 35.705 0
# }