Should co-variates be included in a post-match GLM to estimate treatment effect?

Question

After conducting the following genetic match using the MatchIt package, I would like to conduct a post-match analysis of treatment effect:

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

I can accomplish that with the following:

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)

The weights = weights are weights produced by by the genetic match and appropriately account for the fact that multiple control units are matched to the same treated unit and that each treated unit might have multiple controls.

I have a few questions regarding the post match analysis, feel free to answer as few as you would like:

1. If my outcome is binary, can I simply adjust my model by adding family = "binomial"?

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

2. Do I need to include any additional covariates in the formula outcome ~ hxcopd? My outcome is binomial but I am also interested to hear if the answer to this question depends on if my outcome is either binomial or continuous. Additionally, the explanatory variables are a mix of binomial, categorial, and continuous.

This cross-validated answer suggests that they do not need to be included if balance is obtained by pre-processing - but how can I validate that in my own data quantitatively to say that I am "balanced enough" not to include the covariates. The user then goes on to say doing it anyway is not a bad idea because it is "doubly robust", so does that make this discussion a moot point?

3. If I set replace = FALSE and ties = FALSE so that my genetic matching is essentially nearest neighbor, would I even need the use of glm? Would it be valid just to directly compare the sample populations via hypothesis testing the outcome of interest (such as chi-square)?

Answer 1

The key with binary outcomes is deciding on the causal estimand. First, you need to decide on your effect measure. The most commonly reported effect measures are the risk difference (RD; $P(Y^1=1)-P(Y^0=1)$), the risk ratio (RR; $P(Y^1=1)/P(Y^0=1)$) and the odds ratio (OR; $\frac{P(Y^1=1)}{1-P(Y^1=1)}/\frac{P(Y^0=1)}{1-P(Y^0=1)}$). Second, you need to decide whether you want a conditional effect or a marginal effect. A conditional effect is the effect for a given subclass of units (i.e., defined by a combination of covariate values). A marginal effect is the effect for a population. Conditional effects are useful when deciding whether to treat an individual patient. Marginal effects are useful when deciding whether to implement a rule or policy.

For the RD, the marginal effect is the average of the conditional effects. In this way, the RD is collapsible. The RR and OR are not collapsible, meaning if there is any heterogeneity in individual causal effects, the average of the conditional effects is not the marginal effect (unless the marginal and conditional effects are all zero). This has led to many debates about the best choice of estimand.

Conditional Effects

There will certainly be effect heterogeneity in the conditional RD, which is to say that the RD of the treatment is not the same for all units. For example, if the average RD of for treatment was .4, what would happen to units who had a baseline risk (i.e., under control) of .7? Clearly their risk under treatment cannot be 1.1. If you want to estimate conditional RDs, you have to model the effect heterogeneity, which is sure to be nonlinear. This is true of the RR as well (i.e., how would an average RR of 2 apply to someone with a baseline risk of .55?).

The OR does not face this problem; it is possible for the OR to be the same for all units. No matter what the baseline risk is, the expected risk under treatment is a valid number even if the OR is constant for all units. This is a reason to prefer the OR; you don't need to model effect heterogeneity, and it's possible to provide a single number that represents the conditional effect for all units. This is not possible with the RD or RR. Estimating the conditional OR is also straightforward; it is the coefficient on the treatment in a logistic regression of the outcome on the covariates and treatment, and the Wald or likelihood ratio test for the coefficient is a test of the null hypothesis of no treatment effect.

An often-cited problem with the conditional OR is that its interpretation depends on the variables conditioned on in the model. The results of two studies examining the same treatment and outcome but conditioning on different covariates are not comparable to each other, even in the absence of confounding. In addition, the OR itself is challenging for people to interpret; most people tend to think about risks rather than odds, and people often misinterpret the OR as a RR (which is actually sometimes, but not often, valid).

Marginal Effects

A marginal effect integrates over possible effect heterogeneity in the population and cannot be said to apply to any individual unit when effect heterogeneity is present. That is, an estimated marginal RD of .2 doesn't actually mean the RD is .2 for any individual, which means you can't use a marginal effect to make clinical decisions. Even if there is no effect heterogeneity as measured by the OR, the marginal OR is not equal to the OR for an individual. When there is effect heterogeneity, the marginal effect integrates across the population, but the specific characteristics of the population (which are often unmeasured) determine the marginal effect, so different populations will generally have different marginal effects, and it's rare that a study sample is actually a simple random sample from the population to which the effect is to generalize to. (There are methods of generalizing an effect to a specific population, though, but those are beyond simply estimating an effect from that data you have.)

In a randomized trial, the unadjusted contrast of the outcome under treatment and control is a marginal effect estimate. Because propensity score methods seek to emulate the conditions of a randomized trial (i.e., by removing confounding due to the measured covariates assuming strong ignorability), the unadjusted contrast of the outcome under treatment and control is a marginal effect estimate. Estimating the unadjusted marginal effect in a matched sample is as simple as regressing the outcome on just the treatment using the appropriate link function. (It's possible to estimate a covariate-adjusted marginal effect as well; this is called g-computation and involves average across model-predicted values of the outcome under treatment and control.)

---

Now let's get into the practical aspect of this for you and your analysis. You need to decide what kind of estimand you want. Given that your treatment is COPD, which is not something about which policy can be made, it is likely you are interested in estimating a conditional effect (i.e., the effect of COPD for each unit). If you want a single number that can plausibly represent the conditional effect for all units, you want a conditional OR. If you just want to test the null hypothesis that there is no effect (i.e., what you would do with a $\chi^2$ test), it doesn't matter what method you use, but a covariate-adjusted conditional OR will have the most power to detect the effect.

To get a conditional OR estimate, you need to regress the outcome on the treatment and conditioning covariates using a logistic regression model. We can do this in R with the following code:

fit <- glm(outcome ~ hxcopd + sex + race + age + diabetes + 
              smoke + dyspnea + fnstatus2 + ascites + 
              hxchf + hxmi + hxangina + hypermed + renafail + 
              dialysis + steroid + bleeddis + wtloss, 
           data = df_m2, weights = weights,
           family = quasibinomial(link = "logit"))

lmtest::coeftest(fit, vcov. = sandwich::vcovHC)[1:2,]

The coefficient on treatment is the conditional OR and can be reported as such in the results write-up, but it's important to report the specific covariates you conditioned on to provide context for interpreting the results.

You may ask, if you were going to perform regression on the covariates anyway, why do matching? Matching makes it more likely that a logistic regression model on main effects will be a good fit to your data and prevents minor additional confounding due to unmodelled nonlinearities. That said, some argue that if you are going to perform logistic regression, matching is not necessary as long as work is done to ensure the regression is robust. Frank Harrell's RMS book and course explain this in detail, and his rms package provides tools to do this.

The code you ran estimates a marginal effect (and without family = quasibinomial(link = "logit")), that effect is on the RD scale. While such an analysis would be appropriate for some research questions, it doesn't seem like it would be appropriate for yours.