---
title: "Treatment effects and matching"
date: "2019-01-10"
---
## Treatment effects in observational studies
Most economic data are observational. With unobserved confounders we need instrumental variables, but instruments are hard to find and harder to justify. Under conditional independence — no unobserved confounders given the covariates — we can proceed without them.
Stata's `teffects` suite implements the standard estimators under conditional independence. (The companion `eteffects` handles endogenous treatment via a control-function approach with instruments.) This chapter walks through how `teffects` works.
### Regression adjustment (RA)
RA allows all covariate effects to differ between treatment and control — equivalent to an outcome model with treatment interacted with all covariates.
```{r}
#| label: stata-chunk1
#| engine: 'stata'
#| engine.path: '/usr/local/bin/stata'
#| cache: true
clear
webuse bweightex
teffects ra (bweight prenatal1 mmarried mage fbaby) (mbsmoke)
reg bweight i.mbsmoke##c.(prenatal1 mmarried mage fbaby)
margins r.mbsmoke
```
The data here are `bweightex`, a 60-observation extract with birthweight as the outcome and maternal smoking as the treatment.
`teffects ra` returns the same ATE as `margins r.mbsmoke` after the interacted regression: both give $-389.3099$ grams. The standard errors differ slightly, 37.01 against 36.14, because `teffects` accounts for the estimation of the outcome model differently from `margins`' delta method. The point estimate is the same object computed two ways.
For the ATET (average treatment effect on the treated):
```{r}
#| label: stata-chunk2
#| engine: 'stata'
#| engine.path: '/usr/local/bin/stata'
#| cache: true
clear
webuse bweightex
teffects ra (bweight prenatal1 mmarried mage fbaby) (mbsmoke), atet
reg bweight i.mbsmoke##c.(prenatal1 mmarried mage fbaby)
margins r.mbsmoke, subpop(mbsmoke)
```
The ATET compares potential outcomes for the treated subpopulation, hence the `subpop(mbsmoke)` option in `margins`. Both routes give $-437.4721$, again identical, and again with slightly different standard errors (53.44 against 51.10). Note the ATET is 48 grams larger in magnitude than the ATE of $-389.31$: smoking mothers are drawn disproportionately from the groups where smoking does most damage.
### Inverse Probability Weighting (IPW)
IPW has two steps. First, estimate the treatment model (typically logit) to obtain the propensity score. Second, use inverse-probability weights to compute the outcome difference between treated and control units.
These steps produce consistent estimates of the effect parameters
because the treatment is assumed to be independent of the potential
outcomes after conditioning on the covariates.
We can manually conduct the two steps, but the nice thing about using Stata's teffects is that it takes account of the noise of estimating probability in the first step when calculating standard errors in the second step.
#### example
```{r}
#| label: stata-chunk3
#| engine: 'stata'
#| engine.path: '/usr/local/bin/stata'
#| cache: true
clear
webuse cattaneo2
teffects ipw (bweight) (mbsmoke mmarried c.mage##c.mage fbaby medu, probit)
probit mbsmoke mmarried c.mage##c.mage fbaby medu
predict ps
replace ps = 1/ps if mbsmoke==1
replace ps = 1/(1-ps) if mbsmoke==0
reg bweight mbsmoke [pweight=ps]
```
We switch here to `cattaneo2`, the full 4,642-birth dataset, and a richer propensity model (a probit in marital status, a quadratic in maternal age, first-baby status and education).
`teffects ipw` gives an ATE of $-230.6886$ with a standard error of 25.82. The manual two-step replication matches the point estimate exactly, $-230.6886$; the only difference is the standard error (the manual version does not account for first-stage estimation noise). Note this is well short of the $-389$ found by RA above — but on different data and a different covariate set, so the two are not comparable.
### Doubly robust estimators
RA models the outcome; IPW models the treatment assignment. Doubly-robust estimators model both, requiring only one to be correctly specified.
Stata implements the AIPW combination of Robins and Rotnitzky (1995) and the IPWRA combination of Wooldridge (2010).
The AIPW estimator augments the IPW estimator with a correction term built from the outcome model. Double robustness is best stated as a **consistency** property: AIPW is consistent if *either* the treatment model *or* the outcome model is correctly specified (not necessarily both). It is *not* that the correction term literally "goes to zero" when the treatment model is right and the outcome model is wrong — with a misspecified outcome model the correction term does not generally vanish. Rather, the two pieces are constructed so that the estimator's bias cancels in the limit as long as one of the two models is right.
The IPWRA estimator uses IPW probability weights when performing RA. The weights do not affect the accuracy of the RA estimator if the treatment model is wrong and the outcome model is correct. The weights correct the RA estimator if the treatment model is correct and the outcome model is wrong.
Here we run the two commands, then manually replicate both as a two-step ("plug-in") procedure. `teffects` estimates the treatment and outcome models jointly by GMM/ML, using the joint estimation to get standard errors that account for first-stage estimation noise; the manual two-step version below gets essentially the same point estimate but not the correct standard errors (for the same reason noted above for IPW).
```{r}
#| label: stata-chunk4
#| engine: 'stata'
#| engine.path: '/usr/local/bin/stata'
#| cache: true
clear
webuse cattaneo2
teffects ipwra (bweight mmarried mage prenatal1 fbaby) (mbsmoke mmarried mage prenatal1 fbaby)
teffects aipw (bweight mmarried mage prenatal1 fbaby) (mbsmoke mmarried mage prenatal1 fbaby)
```
**Manual IPWRA:** first estimate the propensity score, then run *weighted* regression adjustment within each treatment arm, using the inverse probability of the *observed* treatment status as weights:
```{r}
#| label: stata-chunk4b
#| engine: 'stata'
#| engine.path: '/usr/local/bin/stata'
#| cache: true
clear
webuse cattaneo2
logit mbsmoke mmarried mage prenatal1 fbaby
predict ps
gen wgt1 = mbsmoke/ps
gen wgt0 = (1-mbsmoke)/(1-ps)
reg bweight mmarried mage prenatal1 fbaby if mbsmoke==1 [pweight=wgt1]
predict mu1_ipwra
reg bweight mmarried mage prenatal1 fbaby if mbsmoke==0 [pweight=wgt0]
predict mu0_ipwra
gen tau_ipwra_i = mu1_ipwra - mu0_ipwra
sum tau_ipwra_i
```
The mean of `tau_ipwra_i` matches `teffects ipwra`'s ATE exactly: $-238.7679$ in both.
**Manual AIPW:** estimate outcome models $\hat\mu_1(X)$, $\hat\mu_0(X)$ (unweighted, one per arm) and the propensity score $\hat p(X)$, then combine via the augmented-IPW formula
$$ \hat\tau_{i}^{aipw} = D_i\frac{Y_i-\hat\mu_1(X_i)}{\hat p(X_i)} + \hat\mu_1(X_i) - (1-D_i)\frac{Y_i-\hat\mu_0(X_i)}{1-\hat p(X_i)} - \hat\mu_0(X_i) $$ {#eq-treatment-matching-1}
```{r}
#| label: stata-chunk4c
#| engine: 'stata'
#| engine.path: '/usr/local/bin/stata'
#| cache: true
clear
webuse cattaneo2
logit mbsmoke mmarried mage prenatal1 fbaby
predict ps
reg bweight mmarried mage prenatal1 fbaby if mbsmoke==1
predict mu1_aipw
reg bweight mmarried mage prenatal1 fbaby if mbsmoke==0
predict mu0_aipw
gen aipw_i = mbsmoke*(bweight-mu1_aipw)/ps + mu1_aipw - (1-mbsmoke)*(bweight-mu0_aipw)/(1-ps) - mu0_aipw
sum aipw_i
```
The mean of `aipw_i` is $-239.0294$, matching `teffects aipw`'s ATE to four
decimals. (One might expect a small gap, since `teffects` estimates the
treatment and outcome models jointly while this is a two-step plug-in; here the
point estimates coincide, and only the standard errors would differ.)
The two doubly-robust estimators land within 0.3 grams of each other,
$-238.77$ for IPWRA and $-239.03$ for AIPW, both about 8 grams away from plain
IPW's $-230.69$.
The unit-level `aipw_i` values are worth a glance before trusting the mean. They
have a standard deviation of 1,618.6 and run from $-29{,}389$ to $+16{,}313$ —
two orders of magnitude wider than the estimate itself. That spread comes from
dividing residuals by $\hat p$ or $1-\hat p$ for units with extreme propensity
scores. The average is well behaved, but it is an average over some very large
individual terms, which is exactly why AIPW can be unstable when overlap is
poor.
## Matching
Matching does not solve the endogeneity problem. It addresses selection on observables only; with unobserved confounders it does not help.
What matching does is balance the covariate distributions across treatment and control groups. Regression does this too, but matching makes the lack of overlap visible — a plain regression can extrapolate far outside the support without warning. Running regression on a matched sample combines the strengths of both approaches.
Below we cover three matching methods in Stata: propensity score matching, nearest-neighbor matching, and coarsened exact matching.
### Propensity score matching
Exact matching on all covariates is ideal but infeasible in high dimensions. Rosenbaum and Rubin (1983) showed that the propensity score provides a scalar sufficient statistic, reducing the matching problem to one dimension.
Stata's `teffects psmatch` and `psmatch2` both estimate the propensity score, match, and compute the treatment effect with correct standard errors in one step.
```{r}
#| label: stata-chunk5
#| engine: 'stata'
#| engine.path: '/usr/local/bin/stata'
#| cache: true
clear
webuse cattaneo2
teffects psmatch (bweight) (mbsmoke mmarried c.mage##c.mage fbaby medu), atet
psmatch2 mbsmoke mmarried c.mage##c.mage fbaby medu, logit ties
reg bweight mbsmoke [aweight=_weight]
```
`teffects psmatch` gives the ATT of smoking on birthweight, $-236.7848$ with a standard error of 26.58. The manual replication uses `psmatch2` with `logit` and `ties` (by default `teffects` includes all tied matches, while `psmatch2` keeps only one), then runs a weighted regression, and gives $-236.7848$ — the same to four decimals. Point estimates match; standard errors should be taken from `teffects`. The nearest-neighbour ATT sits close to the two doubly-robust ATEs above, $-238.8$ and $-239.0$, though it is a different estimand on a different subpopulation.
### Nearest neighbor matching
`teffects nnmatch` uses Mahalanobis distance by default and can also enforce exact matching on discrete covariates. It is nonparametric.
```{r}
#| label: stata-chunk6
#| engine: 'stata'
#| engine.path: '/usr/local/bin/stata'
#| cache: true
clear
webuse cattaneo2
teffects nnmatch (bweight mmarried mage fbaby medu) (mbsmoke) , atet
teffects psmatch (bweight) (mbsmoke mmarried mage fbaby medu), atet
```
Note the specification difference: `nnmatch` lists covariates in the outcome parentheses, `psmatch` lists them in the treatment parentheses.
### Coarsened Exact Matching (CEM)
CEM is available as a separate Stata package.
```{r}
#| label: stata-chunk7
#| engine: 'stata'
#| engine.path: '/usr/local/bin/stata'
#| cache: true
clear
webuse cattaneo2
cem mmarried mage fbaby medu, treatment(mbsmoke)
reg bweight mbsmoke [aweight=cem_weights]
* Note: cem creates the variable "cem_weights" (plural). We spell it out
* in full here rather than relying on Stata's automatic variable-name
* abbreviation matching, to keep the reference unambiguous.
```
CEM creates exact-match strata on coarsened covariates and generates weights for the final regression. The standard errors do not account for the matching step.
## Modern matching in R: balance diagnostics and entropy weighting
The 2019 vintage of this chapter used Stata's `teffects` throughout. R's
`MatchIt` and `WeightIt` packages have since become the standard for matching
workflows, and `cobalt` provides publication-quality balance tables and plots
that are now expected in applied papers.
### Balance checking with cobalt
After any matching or weighting step, always inspect covariate balance.
`cobalt`'s `bal.tab()` and `love.plot()` are the go-to tools:
```{r}
#| echo: true
#| message: false
#| warning: false
#| cache: true
library(MatchIt)
library(cobalt)
# Propensity score matching on the lalonde data
data("lalonde", package = "MatchIt")
m_out <- matchit(treat ~ age + educ + race + married + nodegree +
re74 + re75,
data = lalonde, method = "nearest", distance = "glm",
ratio = 1)
# Balance table: standardized mean differences before and after matching
bal.tab(m_out, stats = c("mean.diffs", "variance.ratios"), thresholds = c(m = 0.1))
# Love plot: visualize balance improvement
love.plot(m_out, stats = "mean.diffs", threshold = 0.1,
abs = TRUE, var.order = "unadjusted")
```
The love plot shows each covariate's standardized mean difference (SMD)
before matching (open circles) and after (filled circles). The vertical line
at SMD = 0.1 is the conventional threshold for "balanced." If post-matching
points cross the line, the matching is not improving balance and a different
specification is needed.
That is exactly what happens here, which is why it is worth running the
diagnostic rather than assuming. 1:1 nearest-neighbour matching on the
propensity score leaves four of the nine covariates *above* the 0.1 threshold,
`race` worst of all (an adjusted SMD of 0.37 for the black indicator), while
discarding 244 of the 429 controls. Matching on a propensity score balances the
score, not necessarily the covariates that went into it.
### Entropy balancing
Propensity score matching discards unmatched units and can reduce effective
sample size substantially. **Entropy balancing** (Hainmueller 2012) keeps all
units and reweights the control group so that weighted covariate moments
exactly match the treatment group — without discarding anyone.
```{r}
#| echo: true
#| message: false
#| warning: false
#| cache: true
library(WeightIt)
library(cobalt)
# Entropy balancing: exact balance on means, variances, and covariances
w_out <- weightit(treat ~ age + educ + race + married + nodegree +
re74 + re75,
data = lalonde, method = "ebal", estimand = "ATT",
moments = 1) # moments=1: balance means only
# moments=2: balance means + variances
summary(w_out)
# Check balance
bal.tab(w_out, stats = "mean.diffs", thresholds = c(m = 0.1))
love.plot(w_out, threshold = 0.1, abs = TRUE)
# Weighted outcome regression
library(lmtest); library(sandwich)
m_eb <- lm(re78 ~ treat + age + educ + race + married + nodegree +
re74 + re75,
data = lalonde, weights = w_out$weights)
coeftest(m_eb, vcov = vcovHC(m_eb, type = "HC3"))["treat", ]
```
Entropy balancing delivers on its promise here: the largest absolute
standardized mean difference after weighting is **exactly zero**, since the
weights are chosen to satisfy the mean-balance constraints as equalities. All
429 controls are retained.
| Method | Pros | Cons |
|---|---|---|
| PSM (nearest neighbor) | Intuitive; respects overlap | Discards units; balance not guaranteed |
| CEM | Exact balance in bins | Requires discretizing continuous vars |
| Entropy balancing | Exact mean balance; keeps all units | Weights can be extreme; no variance balance by default |
| IPW (logistic) | Simple; relies only on the treatment model | Not doubly robust on its own; extreme weights if overlap is poor |
| AIPW / IPWRA | Doubly robust: consistent if *either* the treatment or outcome model is correct | Needs both models specified; more moving parts |
For most applications, entropy balancing or IPW with trimmed weights (via
`WeightIt`'s `trim()`) is preferred over PSM: better effective sample size,
guaranteed balance, and easy integration with doubly-robust estimators.
See also the [weighting chapter](weighting-part2.qmd) for IPW estimators and
the [OLS-ATE chapter](ols-ate.qmd) for regression-based approaches.
---
<!-- see-also-footer -->
*Systematic treatment: [R](https://xiangao.github.io/causal_econometrics_guide/matching.html) · [Julia](https://xiangao.github.io/causal_econometrics_julia/matching.html).*
