Adjusting the model by propensity scores after propensity score matching

Question

I want to control multiple confounders in my data, and I have noticed that including the propensity scores as a variable in the model gives good performance after propensity score matching. I know that including the covariates in the model after propensity scores matching is an established method called doubly robust method, but is it acceptable to use propensity scores instead?

Answer 1

This is definitely possible, and its effectiveness is described in Austin (2017). In general, there is little utility to doing this. If the matching was good, the propensity score distributions should be almost identical between the matched treated and control groups, in which case further adjusting for the propensity score does nothing to reduce bias and will do little to increase precision. If the matching was't good, propensity score adjustment can reduce bias, but it would be better to just try different matching methods until you find one with good balance. If this is not possible, then you will have to make a tradeoff between bias and extrapolation (by using regression on the covariates) or generalizability (by using a method that changes the estimand, like cardinality matching or caliper matching).

If you do decide to do propensity score adjustment after matching, here is what you should do. First, generate splines or polynomials of the logit of the propensity score (to allow the outcome to be flexibly modeled). Then, perform g-computation using the adjusted propensity score. You can use a marginal effects procedure to do this, as demonstrated using the R code below.

library(MatchIt)
data("lalonde")

m.out <- matchit(treat ~ age + educ + race + married +
                   nodegree + re74 + re75,
                 data = lalonde, method = "full")

m.data <- match.data(m.out, distance = "ps")

fit <- lm(re78 ~ treat * splines::ns(qlogis(ps), 5),
          data = m.data, weights = weights)

library(marginaleffects)
g.comp <- comparisons(fit, newdata = subset(m.data, treat == 1),
                      variables = "treat", vcov = ~subclass)
summary(g.comp)
#> Average contrasts 
#>   treat Effect Std. Error z value  Pr(>|z|) 2.5 % 97.5 %
#> 1 1 - 0   1964        719   2.731 0.0063115 554.5   3373
#> 
#> Model type:  lm 
#> Prediction type:  response

^{Created on 2022-06-27 by the reprex package (v2.0.1)}

First, we do full matching on the propensity score using MatchIt::matchit(). Then we extract the matched dataset using match.data(), renaming the propensity score (distance) variable to "ps". Then we fit a model for the outcome with the treatment, a 5 df natural cubic spline, and their interaction as predictors in the matched dataset, including the matching weights. Next, we use marginaleffects::comparisons() to compute the average marginal effects in the treated group (i.e., the ATT) by supplying the model we just fit, the treated units in the matched dataset, and the matched pair membership variable (subclass) for cluster-robust inference to account for the matching. We find an ATT estimate of 1964 with a standard error of 719.

If you were to do covariate adjustment instead, you would replace the linear regression with one that includes the covariates instead of the propensity score, e.g.,

fit <- lm(re78 ~ treat * (age + educ + race + married + nodegree + re74 + re75),
          data = m.data, weights = weights)

A few things are critical here but are often missed. First, it is critical you interact the treatment with the propensity score or covariates in the outcome model. This is necessary to ensure the estimate corresponds to the ATT. It is critical that the matching weights are included in the outcome model. This is necessary to ensure the model is fit in the matched sample. It is critical that cluster-robust standard errors are used. Cluster-robust standard errors are necessary to ensure valid inference after matching.

This procedure would be identical if you had a binary outcome and you wanted the marginal risk difference except that you would replace lm(.) with glm(., family = quasibinomial).