Combine multiple imputation with robust standard errors?

[LukasWallrich]

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.

[bwiernik]

Got a link to that blog post?

[LukasWallrich]

Sorry - https://www.jepusto.com/mi-with-clubsandwich/

[strengejacke]

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}