Why does propensity score matching fail to estimate the true causal effect when OLS works?

Question

Consider the following model (DAG), where D is the treatment (exposure) and Y1 is the outcome. To estimate the causal effect of D on Y1, we can simply condition on V and Y0 (ignoring E).

My issue is that a simple OLS model recover the true model, but when I try to use matching it fails. Why?

Below I simulate the model.

library(tidyverse)
library("MatchIt")
library("optmatch")     
library("Matching")

set.seed(123)
N = 5000

E0 = rnorm(N, 0, 1)
V = rnorm(N, 0, 2)

Y0 = E0*3 + rnorm(N, 0, 2)
D = Y0*5 + (-2*E0) + V*3 + rnorm(N, 0, 5)

# create binary treatment #
Dbin = rbinom(N, 1, plogis(D))

# outcome
Y1 = Y0*6 + (-10*Dbin) + V*2 + rnorm(N, 0, 10)

# put into a dataframe
df = data.frame(E0, Y0, V, D, Dbin, Y1)

# correct model #
# D = -10, causal effect should be approx. -10

Using a simple OLS regression, we recover the true effect of D, which I set to -10

lm(Y1 ~ Dbin + Y0 + V, data = df) %>% summary() 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.26253    0.26069   1.007    0.314    
Dbin        -10.38234    0.43763 -23.724   <2e-16 ***

However, now when I try a Propensity Score Matching (or even other types of matching), I never get to retrieve the true effect.

Here I use the variables Y0 and V because they should be sufficient (using E does not change the estimation).

# Estimate the propensity model
glm1 <- glm(Dbin ~ Y0 + V, family=binomial, data=df)

X  <- glm1$fitted
Y  <- df$Y1
Tr  <- df$Dbin

rr  <- Match(Y=Y, Tr=Tr, X=X, M=1, estimand = "ATT");

I get an ATT of -7, and an ATE of -5!

When I try other matching methods, I also get weird results.

m1 = matchit(Dbin ~ Y0 + V + E0, df, distance = 'mahalanobis', method = 'full')
mat1 = match.data(m1)
lm(Y1 ~ Dbin, mat1, weights = weights) %>% summary()

What am I doing wrong?

Answer 1

As @CloseToC mentioned in the comments, this is because you have a nearly pathological data scenario here. There are a few things that make this scenario "unfair" to matching (i.e., not suitable for matching but well suited for regression). The greatest is that there is essentially no overlap in the propensity score distribution. This is a plot of the true propensity scores between the treatment groups:

There is no way matching, which relies on units of different groups having similar propensity scores, could ever hope to estimate the effect correctly in any population. Using the estimated propensity scores, the story is not much better, and the propensity scores are estimated essentially incorrectly because the distribution is not as extreme as that of the true propensity scores:

There is still a significant lack of overlap. When you perform standard matching (with replacement, as in Matching), almost every treated unit is matched to the very few control units with estimated propensity scores close 1. Indeed, the effective sample size (ESS) of the control group after matching for the ATT is less than 4 (out of an original control sample of 2477). If we look at covariate balance after matching for the ATT, we see significant imbalance remaining in the covariates:

> cobalt::bal.tab(df[c("E0", "V", "Y0")], treat = Tr, 
                  weights = cobalt::get.w(rr), un = T, 
                  method = "m")
Balance Measures
      Type Diff.Un Diff.Adj
E0 Contin.  1.3866  -0.0876
V  Contin.  0.5777   0.1847
Y0 Contin.  2.0520   0.3426

Sample sizes
                     Control Treated
All                  2477.      2523
Matched (ESS)           3.45    2523
Matched (Unweighted)  928.      2523
Unmatched            1549.         0

Let's use a matching method that is actually equipped to deal with poor overlap: matching with a caliper. When we set a very small caliper, so only treated units with control units within their caliper widths are matched and the rest are dropped, we actually get good balance:

> m <- matchit(Dbin ~ Y0 + V + E0, data = df, caliper = .01)
> cobalt::bal.tab(m)
Call
 matchit(formula = Dbin ~ Y0 + V + E0, data = df, caliper = 0.01)

Balance Measures
             Type Diff.Adj
distance Distance   0.0071
Y0        Contin.   0.0175
V         Contin.  -0.0306
E0        Contin.  -0.0048

Sample sizes
          Control Treated
All          2477    2523
Matched       437     437
Unmatched    2040    2086

And when estimating the treatment effect in this sample, we actually get the right answer (because in this case the treatment effect is constant):

> summary(lm(Y1 ~ Dbin, data = match.data(m)))

Call:
lm(formula = Y1 ~ Dbin, data = match.data(m))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.2705     0.5619   0.481     0.63    
Dbin         -9.0438     0.7946 -11.381   <2e-16 ***
---

See also this question, which is almost identical and for which the solution is the same. I will reiterate what I said in that post: treating propensity score matching as a blunt instrument you can just apply without looking at the data is the wrong way to use it. You need to tailor the matching method to the data scenario at hand. By looking at the data, we could see that we were in a low-overlap scenario, so methods that explicitly deal with low overlap should be used. In this case, using caliper matching was successful, but other matching methods like coarsened exact matching (CEM) and cardinality matching, both available in MatchIt, are able to return the correct answer. Methods not well suited for low overlap, like nearest neighbor matching without a caliper or full matching, will not be successful, as you witnessed.