Combine multiple imputation with robust standard errors?
LukasWallrich opened this issue · 4 comments
I am looking for a way to pool multiple imputations (of metafor
models) while using cluster-robust standard errors. Given that parameters
does a great job at each of the two parts, I was hoping to be able to combine the functions - but that does not work.
The procedure for that does not seem to be well established - so far, I only found a [blog post by James E. Pustejovsky] with code - so you might not want to include such a functionality before it is better validated ... but I still wanted to raise the suggestion.
Got a link to that blog post?
Can you elaborate why it doesn't work? I just checked, and it seems fine:
library(mlmRev)
library(mice)
library(dplyr)
library(clubSandwich)
library(parameters)
data(bdf)
bdf <- bdf %>%
select(schoolNR, IQ.verb, IQ.perf, sex, ses, langPRET, aritPRET, aritPOST) %>%
mutate(
schoolNR = factor(schoolNR),
sex = as.numeric(sex)
) %>%
filter(as.numeric(schoolNR) <= 30) %>%
droplevels()
bdf_missing <-
bdf %>%
select(-schoolNR) %>%
ampute(run = TRUE)
bdf_missing <-
bdf_missing$amp %>%
mutate(schoolNR = bdf$schoolNR)
Impute_bdf <- mice(bdf_missing, m = 10, meth = "norm.nob", seed = 24)
models <- lapply(1:5, function(i) {
lm(aritPOST ~ aritPRET + langPRET + sex + ses, data = mice::complete(Impute_bdf, action = i))
})
# not robust
pool_parameters(models)
#> # Fixed Effects
#>
#> Parameter | Coefficient | SE | 95% CI | Statistic | df | p
#> -------------------------------------------------------------------------------
#> (Intercept) | -2.99 | 1.12 | [-5.20, -0.78] | -2.67 | 248.73 | 0.008
#> aritPRET | 1.02 | 0.07 | [ 0.88, 1.17] | 14.19 | 98.90 | < .001
#> langPRET | 0.27 | 0.04 | [ 0.19, 0.35] | 7.01 | 67.77 | < .001
#> sex | -0.84 | 0.44 | [-1.73, 0.05] | -1.89 | 55.81 | 0.064
#> ses | 0.02 | 0.02 | [-0.02, 0.07] | 1.19 | 67.91 | 0.238
#>
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#> using a Wald distribution approximation.
# robust
pool_parameters(models, vcov = "CR2", vcov_args = list(cluster = bdf$schoolNR))
#> # Fixed Effects
#>
#> Parameter | Coefficient | SE | 95% CI | Statistic | df | p
#> -------------------------------------------------------------------------------
#> (Intercept) | -2.99 | 1.36 | [-5.66, -0.32] | -2.21 | 340.15 | 0.028
#> aritPRET | 1.02 | 0.08 | [ 0.87, 1.18] | 12.91 | 131.93 | < .001
#> langPRET | 0.27 | 0.04 | [ 0.20, 0.34] | 7.57 | 51.82 | < .001
#> sex | -0.84 | 0.50 | [-1.83, 0.15] | -1.69 | 82.89 | 0.096
#> ses | 0.02 | 0.03 | [-0.03, 0.08] | 0.82 | 203.05 | 0.413
#>
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#> using a Wald distribution approximation.
Created on 2024-02-28 with reprex v2.1.0
Closing for now. If you still see problems, please re-open or file a new issue.