So?

OK - so maybe that isn’t really all that interesting. Why am I even talking about this? Well, in the actual study, we have a fair amount of missing data. In some cases we don’t have an EK measure, and in others we don’t have an outcome measure. And since the missingness is on the order of 15% to 20%, we decided to use multiple imputation. We used the mice package in R to impute the data, and we pooled the model estimates from the completed data sets to get our final estimates. mice is a fantastic package, but one thing that is does not supply is some sort of pooled variance-covariance matrix. Looking for a relatively quick solution, I decided to use bootstrap methods to estimate the confidence intervals.

(“Relatively quick” is itself a relative term, since bootstrapping and imputing together is not exactly a quick process - maybe something to work on. I was also not fitting standard linear models but mixed effect models. Needless to say, it took a bit of computing time to get my estimates.)

Seeking credit (and maybe some sympathy) for all of my hard work, I mentioned this laborious process to my collaborator. She told me that you can easily estimate the group specific effects merely by changing the reference group and refitting the model. I could see right away that the point estimates would be fine, but surely the standard errors would not be estimated correctly? Of course, a few simulations ensued.

First, I just changed the reference group so that \(\beta_2\) would be measuring the relationship of EK and academic performance for immigrant kids, and \(\beta_2 + \beta_3\) would represent the relationship for the non-immigrant kids. Here are the levels before the change:

head(dT$fI)
## [1] Imm     Imm     not Imm Imm     not Imm not Imm
## Levels: not Imm Imm

And after:

dT$fI <- relevel(dT$fI, ref="Imm")
head(dT$fI)
## [1] Imm     Imm     not Imm Imm     not Imm not Imm
## Levels: Imm not Imm

And the model:

lm2 <- lm(T ~ fI*EK, data = dT)
tidy(lm2)
##           term  estimate std.error  statistic       p.value
## 1  (Intercept) 11.778770 0.2086526  56.451588 8.367177e-143
## 2    fInot Imm -1.616929 0.2641919  -6.120281  3.661090e-09
## 3           EK  2.065308 0.1045418  19.755813  1.112574e-52
## 4 fInot Imm:EK -1.603680 0.1396076 -11.487049  9.808529e-25

The estimate for this new \(\beta_2\) is 2.07 (se=0.105), pretty much aligned with our estimate that required a little more work. While this is not a proof by any means, I did do variations on this simulation (adding other covariates, changing the strength of association, changing sample size, changing variation, etc.) and both approaches seem to be equivalent. I even created 10000 samples to see if the coverage rates of the 95% confidence intervals were correct. They were. My collaborator was right. And I felt a little embarrassed, because it seems like something I should have known.

But …

Would this still work with missing data? Surely, things would go awry in the pooling process. So, I did one last simulation, generating the same data, but then added missingness. I imputed the missing data, fit the models, and pooled the results (including pooled 95% confidence intervals). And then I looked at the coverage rates.

First I added some missingness into the data

defM <- defMiss(varname = "EK", formula = "0.05 + 0.10*I", 
                logit.link = FALSE)
defM <- defMiss(defM, varname = "T", formula = "0.05 + 0.05*I", 
                logit.link = FALSE)
defM
##    varname       formula logit.link baseline monotonic
## 1:      EK 0.05 + 0.10*I      FALSE    FALSE     FALSE
## 2:       T 0.05 + 0.05*I      FALSE    FALSE     FALSE

And then I generated 500 data sets, imputed the data, and fit the models. Each iteration, I stored the final model results for both models (in one where the reference is non-immigrant and the the other where the reference group is immigrant).

library(mice)

nonRes <- list()
immRes <- list()

set.seed(3298348)

for (i in 1:500) {
  
  dT <- genData(250, def)
  dT <- genFactor(dT, "I", labels = c("non Imm", "Imm"), prefix = "non")
  dT$immI <- relevel(dT$nonI, ref = "Imm")
  
  # generate a missing data matrix
  
  missMat <- genMiss(dtName = dT, missDefs = defM, idvars = "id")
  
  # create obseverd data set
  
  dtObs <- genObs(dT, missMat, idvars = "id")
  dtObs <- dtObs[, .(I, EK, nonI, immI, T)]
  
  # impute the missing data (create 20 data sets for each iteration)

  dtImp <- mice(data = dtObs, method = 'cart', m = 20, printFlag = FALSE)
  
  # non-immgrant is the reference group

  estImp <- with(dtImp, lm(T ~ nonI*EK))
  lm1 <- summary(pool(estImp), conf.int = TRUE)
  dt1 <- as.data.table(lm1)
  dt1[, term := rownames(lm1)]
  setnames(dt1, c("2.5 %", "97.5 %"), c("conf.low", "conf.high"))
  dt1[, iter := i]
  nonRes[[i]] <- dt1
  
  # immgrant is the reference group

  estImp <- with(dtImp, lm(T ~ immI*EK))
  lm2 <- summary(pool(estImp), conf.int = TRUE)
  dt2 <- as.data.table(lm2)
  dt2[, term := rownames(lm2)]
  setnames(dt2, c("2.5 %", "97.5 %"), c("conf.low", "conf.high"))
  dt2[, iter := i]
  immRes[[i]] <- dt2
}

nonRes <- rbindlist(nonRes)
immRes <- rbindlist(immRes)

The proportion of confidence intervals that contain the true values is pretty close to 95% for both estimates:

mean(nonRes[term == "EK", conf.low < 0.5 & conf.high > 0.5])
## [1] 0.958
mean(immRes[term == "EK", conf.low < 2.0 & conf.high > 2.0])
## [1] 0.948

And the estimates of the mean and standard deviations are also pretty good:

nonRes[term == "EK", .(mean = round(mean(estimate),3), 
                       obs.SD = round(sd(estimate),3), 
                       avgEst.SD = round(sqrt(mean(std.error^2)),3))]
##     mean obs.SD avgEst.SD
## 1: 0.503  0.086     0.088
immRes[term == "EK", .(mean = round(mean(estimate),3), 
                       obs.SD = round(sd(estimate),3), 
                       avgEst.SD = round(sqrt(mean(std.error^2)),3))]
##     mean obs.SD avgEst.SD
## 1: 1.952  0.117     0.124

Because I like to include at least one visual in a post, here is a plot of the 95% confidence intervals, with the CIs not covering the true values colored blue:

The re-reference approach seems to work quite well (in the context of this simulation, at least). My guess is the hours of bootstrapping may have been unnecessary, though I haven’t fully tested all of this out in the context of clustered data. My guess is it will turn out OK in that case as well.

Appendix: ggplot code

nonEK <- nonRes[term == "EK", .(iter, ref = "Non-immigrant",
  estimate, conf.low, conf.high,
  cover = (conf.low < 0.5 & conf.high > 0.5))]

immEK <- immRes[term == "EK", .(iter, ref = "Immigrant",
  estimate, conf.low, conf.high,
  cover = (conf.low < 2 & conf.high > 2))]

EK <- rbindlist(list(nonEK, immEK))

vline <- data.table(xint = c(.5, 2), 
                    ref = c("Non-immigrant", "Immigrant"))

ggplot(data = EK, aes(x = conf.low, xend = conf.high, y = iter, yend = iter)) +
  geom_segment(aes(color = cover)) +
  geom_vline(data=vline, aes(xintercept=xint), lty = 3) +
  facet_grid(.~ ref, scales = "free") +
  theme(panel.grid = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        legend.position = "none") +
  scale_color_manual(values = c("#5c81ba","grey75")) +
  scale_x_continuous(expand = c(.1, 0), name = "95% CI")