I´m trying to get up to speed with multiple imputation and matching. I´ve been struggling with finding the correct method of pooling the matched imputed data. In the end I ended up doing a new analysis of the Lalonde data (see below- copy paste, takes approx 60sec to run).
As far as I understand there are two general approaches and detailed by Mitra et al (also see this thread). In the analysis below I did these approaches, but also one where I clustered according to each individual and I obtain the exact same results as the averaging by Mitra (although easier to implement). They seem to fit reasonably well with the original dataset (could probably be improved by increasing the imputations, growing longer trees and increasing depth). Now the big question is which method I should report?
data("lalonde")
lalonde$age[sample(1:614, 50)] <- NA
library(cobalt)
library(mice)
library(tidyverse)
library(twang)
s <- mice(lalonde, m = 3)
s_long <- mice::complete(s, "long")
s_list <- s_long %>% group_by(.imp) %>% nest()
s_list$twang <- map(
s_list$data,
~ ps(
treat ~ age + educ + black + hispan + married + nodegree + re74 + re75,
data = as.data.frame(.),
n.trees = 500,
interaction.depth = 1,
shrinkage = 0.01,
perm.test.iters = 0,
stop.method = c("es.mean", "ks.max"),
estimand = "ATT",
verbose = TRUE
)
)
s_list$w <-
map(s_list$twang,
~ get.weights(
ps1 = .,
stop.method = "es.mean",
estimand = "ATT"
))
s_list$design <- map2(s_list$data,
s_list$w,
~ svydesign(
ids = ~ 1,
weights = ~ .y,
data = .x
))
s_list$glm2 <- map(s_list$design, ~
svyglm(re78 ~ treat, design = .x))
library(broom)
map_df(s_list$glm2, ~ tidy(.), .id = "N")
library(mitools)
summary(MIcombine(s_list$glm2)) #According to Rubins rule
s_long_data <- unnest(s_list, ... = data, w)
glm_id_design <- svydesign(ids = ~ .id,
weights = ~ w,
data = s_long_data)
glm_id <- svyglm(re78 ~ treat,
design = glm_id_design)
##According to Mitra
ave_w <- s_long_data %>% group_by(.id) %>% summarise(w2 = mean(w))
s_mitra <- s_long_data %>% filter(.imp == 1) %>% left_join(ave_w)
glm_mitra_design <- svydesign(ids = ~ 1,
weights = ~ w2,
data = s_mitra)
glm_mitra <- svyglm(re78 ~ treat,
design = glm_mitra_design)
##### All
(glm_all <-
list(
mitra = tidy(glm_mitra),
svy_cluster = tidy(glm_id),
robin = summary(MIcombine(s_list$glm2))
))
##### Compared to
data(lalonde)
all_lalonde <- ps(
treat ~ age + educ + black + hispan + married + nodegree + re74 + re75,
data = lalonde,
n.trees = 500,
interaction.depth = 1,
shrinkage = 0.01,
perm.test.iters = 0,
stop.method = c("es.mean", "ks.max"),
estimand = "ATT",
verbose = TRUE
)
all_svy <- svydesign(
ids = ~ 1,
weights = get.weights(all_lalonde, "es.mean", estimand = "ATT"),
data = lalonde
)
svy_all <- svyglm(re78 ~ treat,
design = all_svy)
(glm_all <- list(
mitra = tidy(glm_mitra),
svy_cluster = tidy(glm_id),
robin = summary(MIcombine(s_list$glm2)),
all_lalonde = tidy(svy_all)
))