Assessing for causality after genetic matching - how to use weights

Question

I am conducting an analysis of the effect of COPD on particular outcomes after surgery. I have found that utilizing the matchit package with the genetic method produces the best balance:

df_match <- MatchIt::matchit(hxcopd ~ sex + race + age + diabetes + smoke + dyspnea + fnstatus2 + ascites + hxchf + hxmi + hxangina + hypermed + renafail + dialysis + steroid + bleeddis + wtloss, data = df_m, method = "genetic", pop.size = 1000)

summary(df_match)

Summary of balance for matched data:
                                        Means Treated Means Control SD Control Mean Diff eQQ Med eQQ Mean eQQ Max
distance                                       0.1637        0.1607     0.1847    0.0030  0.0559   0.1136  0.3928
sexFALSE                                       0.2409        0.2409     0.4280    0.0000  0.0000   0.1168  1.0000
sexTRUE                                        0.7591        0.7591     0.4280    0.0000  0.0000   0.1168  1.0000
raceAsian                                      0.0073        0.0073     0.0852    0.0000  0.0000   0.0000  0.0000
raceBlack                                      0.0876        0.0949     0.2933   -0.0073  0.0000   0.0657  1.0000
raceNative Hawaiian or Pacific islander        0.0000        0.0000     0.0000    0.0000  0.0000   0.0000  0.0000
raceWhite                                      0.8905        0.8905     0.3125    0.0000  0.0000   0.0803  1.0000
age                                           68.2555       68.1509    12.7220    0.1046  1.0000   1.2409  5.0000
diabetesTRUE                                   0.0949        0.0876     0.2829    0.0073  0.0000   0.0730  1.0000
smokeTRUE                                      0.3869        0.3869     0.4874    0.0000  0.0000   0.2044  1.0000
dyspneaTRUE                                    0.3212        0.3212     0.4673    0.0000  0.0000   0.2701  1.0000
fnstatus2Partially dependent                   0.0146        0.0146     0.1200    0.0000  0.0000   0.0073  1.0000
fnstatus2Totally dependent                     0.0000        0.0000     0.0000    0.0000  0.0000   0.0000  0.0000
ascitesTRUE                                    0.0146        0.0146     0.1200    0.0000  0.0000   0.0073  1.0000
hxchfTRUE                                      0.0073        0.0073     0.0852    0.0000  0.0000   0.0000  0.0000
hxmiTRUE                                       0.0000        0.0000     0.0000    0.0000  0.0000   0.0000  0.0000
hxanginaTRUE                                   0.0219        0.0219     0.1465    0.0000  0.0000   0.0146  1.0000
hypermedTRUE                                   0.5766        0.5766     0.4945    0.0000  0.0000   0.0292  1.0000
renafailTRUE                                   0.0073        0.0073     0.0852    0.0000  0.0000   0.0000  0.0000
dialysisTRUE                                   0.0146        0.0073     0.0852    0.0073  0.0000   0.0073  1.0000
steroidTRUE                                    0.1022        0.1022     0.3031    0.0000  0.0000   0.0803  1.0000
bleeddisTRUE                                   0.0657        0.0657     0.2479    0.0000  0.0000   0.0511  1.0000
wtlossTRUE                                     0.0073        0.0073     0.0852    0.0000  0.0000   0.0000  0.0000

Percent Balance Improvement:
                                        Mean Diff. eQQ Med  eQQ Mean  eQQ Max
distance                                   97.8044 19.5737   15.4411  18.9160
sexFALSE                                  100.0000  0.0000 -100.0000   0.0000
sexTRUE                                   100.0000  0.0000 -100.0000   0.0000
raceAsian                                 100.0000  0.0000  100.0000 100.0000
raceBlack                                  52.3051  0.0000 -350.0000   0.0000
raceNative Hawaiian or Pacific islander   100.0000  0.0000  100.0000 100.0000
raceWhite                                 100.0000  0.0000  -83.3333   0.0000
age                                        99.0985 91.6667   89.3149  78.2609
diabetesTRUE                               71.5156  0.0000 -233.3333   0.0000
smokeTRUE                                 100.0000  0.0000    0.0000   0.0000
dyspneaTRUE                               100.0000  0.0000    7.5000   0.0000
fnstatus2Partially dependent              100.0000  0.0000    0.0000   0.0000
fnstatus2Totally dependent                100.0000  0.0000  100.0000 100.0000
ascitesTRUE                               100.0000  0.0000    0.0000   0.0000
hxchfTRUE                                 100.0000  0.0000    0.0000   0.0000
hxmiTRUE                                  100.0000  0.0000  100.0000 100.0000
hxanginaTRUE                              100.0000  0.0000    0.0000   0.0000
hypermedTRUE                              100.0000  0.0000   87.0968   0.0000
renafailTRUE                              100.0000  0.0000    0.0000   0.0000
dialysisTRUE                               17.1281  0.0000    0.0000   0.0000
steroidTRUE                               100.0000  0.0000    0.0000   0.0000
bleeddisTRUE                              100.0000  0.0000  -40.0000   0.0000
wtlossTRUE                                100.0000  0.0000    0.0000   0.0000

Sample sizes:
          Control Treated
All          4490     137
Matched       611     137
Unmatched    3879       0
Discarded       0       0

This seems to be the best balance I can get. I am unsure of the best way to conduct post-matching analysis after a genetic match. How can I best assess the causality of COPD (hxcopd)? I am particularly confused because after a nearest neighbor match, I can just use the resulting dataset to perform whatever analyses I would normally go on to perform.

But the genetic algorithm has balanced the covariates with weights; for example:

If I just check the mean of diabetes in the resulting dataset, it is not concordant with the summary:

df_m2 <- MatchIt::match.data(df_match)

tapply(df_m2$diabetes, df_m2$hxcopd, mean)
     FALSE       TRUE 
0.02291326 0.09489051 

But when multiplied by the weights, it is:

tapply(df_m2$diabetes*df_m2$weights, df_m2$hxcopd, mean)
     FALSE       TRUE 
0.08029197 0.09489051 

I'm therefore not even really sure how to assess the balance or how to go on and perform additional analyses with these weights. For example, can I just directly compare the rates of some outcome (such as infection) between these two groups after matching?

Answer 1

That's some amazing balance! There are a few things you should know about genetic matching with MatchIt. These are due to the fact that MatchIt calls the function GenMatch in the Matching package, which has a different syntax from matchit().

First, by default, it performs matching with replacement, which is not true of nearest neighbor matching. To perform matching without replacement, you need to specify replace = FALSE. Second, by default, it performs variable-ratio matching with ties in that if multiple control units are equally close to a treated unit (i.e., tied), they are all matched to that treated unit. You might think it unusual that there would be so many ties, but whether two units are considered to be tied is controlled by the distance.tolerance option in GenMatch(), which is 1e-5 by default, but could be smaller. You can also set ties = FALSE, which, rather than matching every tied control unit to the corresponding treated unit, randomly selects one of the control units to match.

With replace = FALSE and ties = FALSE, genetic matching is nearest neighbor matching (with balance optimization), and you can estimate the treatment effect in the same way. Otherwise, you have to incorporate weights into the effect estimation and use a robust standard error to account for them. It's good practice to do this even with nearest neighbor matching because including weights and using a robust standard error is compatible with all matching methods. Here's how you would do this:

df_m2 <- MatchIt::match.data(df_match)
fit <- glm(outcome ~ hxcopd, data = df_m2, weights = weights)
lmtest::coeftest(fit, vcov. = sandwich::vcovHC)
lmtest::coefci(fit, vcov. = sandwich::vcovHC)

Setting weights = weights causes glm() (or lm(), or coxph(), etc.) to use the weights stored in the match.data() output. These weights appropriately account for the fact that multiple control units are matched to the same treated unit (if replace = TRUE) and that each treated unit might have multiple controls (if ties = TRUE). If your outcome is continuous and you are using a linear model, you can also include covariates in the outcome model. lmtest provides the functions coeftest() and coefci(), which produce effect estimates, standard errors, and confidence intervals that can incorporate a function to estimate robust standard errors. Using sandwich::vcovHC() uses the "HC3" robust standard error, which is robust to heteroscedasticity and appropriate for small (and large) samples.

---

I'm not one of the original authors of MatchIt, but I'm in the process of updating it after several years without updates. Part of those updates include setting replace = FALSE and ties = FALSE by default to be consistent with other matching methods, so in the future, the output of method = "genetic" will be identical in form to that from method = "nearest". In addition, I've written a vignette detailing how to estimate treatment effects after each type of matching for binary, continuous, and survival outcomes. As of now (2020-09-29), these are available on my GitHub and will likely be on CRAN in the next month or so.