2
$\begingroup$

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)
))
$\endgroup$

1 Answer 1

1
$\begingroup$

The following is my opinion, so feel free to look elsewhere. I think it's really good that you tried several methods and got the same result with each; that shows that your result does not depend on your methods choices, which gives some evidence for a robust effect. I would recommend reporting the method that is the most intuitive, the best-justified, or the most familiar to readers in your field. I think it would be a huge strength to let readers know that using other methods yielded the same result; this can go into an Appendix at most or a footnote at least, but I think a sentence or two in your write-up would be sufficient.

This approach is sometimes used when it isn't clear whether to endorse an assumption, so you try methods that either require or don't require the assumption. If the results agree, you can report either method, but often the simpler method allows readers to arrive at their conclusion the most efficiently. For example, if you didn't want to assume homoscedasticity in ANOVA, you could use Welch's ANOVA, but if you got the same results with the regular uncorrected ANOVA, you might report that instead, and also let readers know that the results didn't differ if you were to use the corrected ANOVA.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.