Tutorial 9: Matching

Evaluation of Public Policies

Edoardo Alberto Viganò

In this lab we practise matching methods using the classic LaLonde (1986) job-training dataset.
We estimate propensity scores, perform nearest-neighbour and caliper matching, assess balance
with love plots, and compare estimates across methods.

1 Setup

# Install once (run this, then comment it out)

# [Link](c("tidyverse", "MatchIt", "cobalt", "optmatch", "sandwich", "lmtest"))

library(tidyverse) # data wrangling and plotting

library(MatchIt) # matching (PSM, Mahalanobis, CEM, etc.)
library(cobalt) # balance tables and love plots
library(sandwich) # robust standard errors
library(lmtest) # coeftest() for robust SEs

2 The data: LaLonde (1986)

This dataset comes from Dehejia and Wahba (1999), a subset of LaLonde’s (1986) National
Supported Work Demonstration. The goal: estimate the effect of a job-training programme on
post-programme earnings.
The experimental benchmark – the effect estimated from the original RCT – is approximately
$1,794. We will try to recover this from observational data using matching.

1
2.1 Loading and inspecting the data

lalonde <- read_csv("[Link]") |>

select(-1) # drop the row-number column

glimpse(lalonde)

Rows: 16,177
Columns: 12
$ treat <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
$ age <dbl> 37, 22, 30, 27, 33, 22, 23, 32, 22, 33, 19, 21, 18, 27, 17, 1~
$ educ <dbl> 11, 9, 12, 11, 8, 9, 12, 11, 16, 12, 9, 13, 8, 10, 7, 10, 13,~
$ black <dbl> 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
$ hispan <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ married <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0~
$ nodegree <dbl> 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1~
$ re74 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ re75 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ re78 <dbl> 9930.0459, 3595.8940, 24909.4492, 7506.1460, 289.7899, 4056.4~
$ unem74 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
$ unem75 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~

2.2 Variable descriptions

Variable Description
treat 1 if assigned to job training, 0 if control
age Age in years
educ Years of education
black 1 if Black
hispan 1 if Hispanic
married 1 if married
nodegree 1 if no school degree
re74 Real earnings in 1974 (pre-treatment)
re75 Real earnings in 1975 (pre-treatment)
re78 Real earnings in 1978 (outcome)
unem74 1 if unemployed in 1974
unem75 1 if unemployed in 1975

2
2.3 How many observations per group?

lalonde |> count(treat)

# A tibble: 2 x 2
treat n
<dbl> <int>
1 0 15992
2 1 185

The dataset is heavily imbalanced: many more control units than treated. This is typical of
observational data and is exactly the situation where matching is useful.

3 Naive estimates (no matching)

3.1 Simple difference in means

model_naive <- lm(re78 ~ treat, data = lalonde)

summary(model_naive)

Call:
lm(formula = re78 ~ treat, data = lalonde)

Residuals:
Min 1Q Median 3Q Max
-14856 -9006 1444 10744 53959

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14855.64 76.22 194.90 <2e-16 ***
treat -8506.50 712.77 -11.93 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9639 on 16175 degrees of freedom

Multiple R-squared: 0.008729, Adjusted R-squared: 0.008668
F-statistic: 142.4 on 1 and 16175 DF, p-value: < 2.2e-16

3
The coefficient is negative, suggesting the programme hurts earnings. But we know from the
experiment that it helps. The problem: treated and control individuals differ systematically on
pre-treatment characteristics. This is a classical example of selection bias.

3.2 Covariate-adjusted regression

options(scipen = 999)

model_ols <- lm(re78 ~ treat + age + educ + black + hispan + nodegree +

married + unem74 + unem75 + re74 + re75, data = lalonde)
summary(model_ols)

Call:
lm(formula = re78 ~ treat + age + educ + black + hispan + nodegree +
married + unem74 + unem75 + re74 + re75, data = lalonde)

Residuals:
Min 1Q Median 3Q Max
-25032 -3524 1294 3774 53742

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5762.17979 445.61449 12.931 < 0.0000000000000002 ***
treat 1067.54611 554.05954 1.927 0.054026 .
age -94.54102 6.00028 -15.756 < 0.0000000000000002 ***
educ 175.22548 28.69658 6.106 0.0000000010440 ***
[ reached getOption("[Link]") -- omitted 8 rows ]
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6999 on 16165 degrees of freedom

Multiple R-squared: 0.4777, Adjusted R-squared: 0.4773
F-statistic: 1344 on 11 and 16165 DF, p-value: < 0.00000000000000022

The estimate moves closer to $1,794, but regression assumes a linear functional form and
extrapolates beyond the data. Matching offers a non-parametric alternative.

4
4 Pre-matching balance

Before matching, we should check how different the treated and control groups are on pre-
treatment covariates. If they are very different, a simple regression comparison is suspect.

[Link](treat ~ age + educ + black + hispan + nodegree +

married + unem74 + unem75 + re74 + re75,
data = lalonde,
stats = "[Link]",
thresholds = c(m = 0.1),
abs = TRUE,
[Link] = "unadjusted",
title = "Covariate balance before matching")

Covariate balance before matching

re75

re74

age

black

educ Sample
unem74 Unadjusted

married

unem75

nodegree

hispan

0.0 0.5 1.0 1.5

Absolute Mean Differences
The dashed line at 0.1 is a common threshold for acceptable balance (standardised mean
differences < 0.1). Covariates beyond this line are poorly balanced. Several covariates are
far from zero, which explains why the naive estimate is very different from the true (i.e.,
experimental) estimate.

5
5 Propensity score estimation

The propensity score is the probability of receiving treatment given pre-treatment covariates:
𝑒(𝑋) = 𝑃 (treat = 1 ∣ 𝑋). Matching on this single number is equivalent to matching on all
covariates simultaneously (Rosenbaum and Rubin 1983).
We estimate it with a logit model:

ps_model <- glm(treat ~ age + educ + black + hispan + nodegree +

married + re74 + re75 + unem74 + unem75,
data = lalonde, family = binomial("logit"))
summary(ps_model)

Call:
glm(formula = treat ~ age + educ + black + hispan + nodegree +
married + re74 + re75 + unem74 + unem75, family = binomial("logit"),
data = lalonde)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -6.34059728 0.81249585 -7.804 0.000000000000006 ***
age -0.01789680 0.01063912 -1.682 0.09254 .
educ 0.01950110 0.04858679 0.401 0.68815
black 4.28584459 0.26262006 16.320 < 0.0000000000000002 ***
[ reached getOption("[Link]") -- omitted 7 rows ]
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 2022.14 on 16176 degrees of freedom

Residual deviance: 950.46 on 16166 degrees of freedom
AIC: 972.46

Number of Fisher Scoring iterations: 10

Several covariates are strong predictors of treatment status, confirming the hypothesis that
assignment is not random.

6
5.1 Visualising the propensity score

lalonde <- lalonde |>

mutate(pscore = predict(ps_model, type = "response"))

ggplot(lalonde, aes(x = pscore, fill = factor(treat))) +

geom_density(aes(y = after_stat(scaled)), alpha = 0.2) +
labs(x = "Propensity score", y = "Density", fill = "Treatment") +
theme_minimal()

1.00

0.75

Treatment
Density

0.50 0
1

0.25

0.00
0.0 0.2 0.4 0.6
Propensity score

If the two densities barely overlap, matching will struggle because there is no good coun-
terfactual for many treated units. This is the common support problem. What do you see
here?

6 Matching with MatchIt

The MatchIt workflow has three steps:

1. Match: use matchit() to pair treated and control units

2. Check balance: did matching make the groups comparable?
3. Estimate: run the outcome model on the matched data

7
6.1 Nearest-neighbour propensity score matching

matchit() matches each treated unit to its nearest control unit (without replacement). In
practice, the non-participant(s) 𝑃𝑗̄ closest to our participant 𝑃𝑖̄ will be selected as a match. The
estimator has the option to perform matching with or without replacement. Today, we will
only use matching without replacement. It’s good to know though that this method has the
disadvantage that the final estimate will usually depend on the initial ordering of the treated
observations for which the matches were selected.

m_nn <- matchit(treat ~ age + educ + black + hispan + nodegree +

married + re74 + re75 + unem74 + unem75,
data = lalonde,
method = "nearest",
distance = "glm", # logistic regression for PS
[Link] = "data",
replace = FALSE)

m_nn

A matchit object
- method: 1:1 nearest neighbor matching without replacement
- distance: Propensity score
- estimated with logistic regression
- number of obs.: 16177 (original), 370 (matched)
- target estimand: ATT
- covariates: age, educ, black, hispan, nodegree, married, re74, re75, unem74, unem75

Let’s check the summary:

summary(m_nn)

Call:
matchit(formula = treat ~ age + educ + black + hispan + nodegree +
married + re74 + re75 + unem74 + unem75, data = lalonde,
method = "nearest", distance = "glm", replace = FALSE, [Link] = "data")

Summary of Balance for All Data:

Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance 0.3217 0.0078 1.3899 27.0244 0.5065
age 25.8162 33.2252 -1.0355 0.4196 0.1863

8
educ 10.3459 12.0275 -0.8363 0.4905 0.0908
eCDF Max
distance 0.8340
age 0.3427
educ 0.4123
[ reached getOption("[Link]") -- omitted 8 rows ]

Summary of Balance for Matched Data:

Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance 0.3217 0.3029 0.0831 1.1002 0.0006
age 25.8162 26.2108 -0.0551 0.4386 0.0839
eCDF Max Std. Pair Dist.
distance 0.1081 0.1060
age 0.2000 1.1113
[ reached getOption("[Link]") -- omitted 9 rows ]

Sample Sizes:
Control Treated
All 15992 185
Matched 185 185
Unmatched 15807 0
Discarded 0 0

The summary shows balance statistics before and after matching. Look at the “Std. Mean Diff.”
column: values should be close to zero after matching.

6.2 Love plot: balance after nearest-neighbour matching

[Link](m_nn,
stats = "[Link]",
thresholds = c(m = 0.1),
abs = TRUE,
[Link] = "unadjusted",
title = "Balance: nearest-neighbour PS matching")

9
 Balance: nearest−neighbour PS matching
distance
re75
re74
age
educ Sample
black Unadjusted

unem74 Adjusted

married
unem75
nodegree
hispan

0 1 2 3
Absolute Mean Differences
Compare this to the pre-matching love plot. Are the covariates now within the 0.1 threshold?

6.3 Estimating the treatment effect

We extract the matched data and run our outcome regression on it:

matched_nn <- [Link](m_nn)

# Simple difference in means on matched data

model_nn <- lm(re78 ~ treat, data = matched_nn)
coeftest(model_nn, vcov. = vcovCL(model_nn, cluster = ~ subclass))

t test of coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 4392.88 427.67 10.2716 < 0.00000000000000022 ***
treat 1956.26 673.07 2.9065 0.003877 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

How does the estimate compare to the experimental benchmark of $1,794?

10
7 Caliper matching

Nearest-neighbour matching can produce bad matches if the closest control is still far away.
Caliper matching imposes a maximum distance: if no control unit is within the caliper, the
treated unit is dropped.

m_caliper <- matchit(treat ~ age + educ + black + hispan + nodegree +

married + re74 + re75 + unem74 + unem75,
data = lalonde,
method = "nearest",
distance = "glm",
caliper = 0.1, # 0.1 SD of the logit PS
[Link] = "data",
replace = FALSE)

m_caliper

A matchit object
- method: 1:1 nearest neighbor matching without replacement
- distance: Propensity score [caliper]
- estimated with logistic regression
- caliper: <distance> (0.006)
- number of obs.: 16177 (original), 298 (matched)
- target estimand: ATT
- covariates: age, educ, black, hispan, nodegree, married, re74, re75, unem74, unem75

Some treated units may be dropped. Check how many:

summary(m_caliper)

Call:
matchit(formula = treat ~ age + educ + black + hispan + nodegree +
married + re74 + re75 + unem74 + unem75, data = lalonde,
method = "nearest", distance = "glm", replace = FALSE, [Link] = "data",
caliper = 0.1)

Summary of Balance for All Data:

Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance 0.3217 0.0078 1.3899 27.0244 0.5065

11
age 25.8162 33.2252 -1.0355 0.4196 0.1863
educ 10.3459 12.0275 -0.8363 0.4905 0.0908
eCDF Max
distance 0.8340
age 0.3427
educ 0.4123
[ reached getOption("[Link]") -- omitted 8 rows ]

Summary of Balance for Matched Data:

Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance 0.2703 0.2704 -0.0005 0.9985 0.0001
age 25.9060 25.5503 0.0497 0.6272 0.0478
eCDF Max Std. Pair Dist.
distance 0.0268 0.0046
age 0.1477 0.9915
[ reached getOption("[Link]") -- omitted 9 rows ]

Sample Sizes:
Control Treated
All 15992 185
Matched 149 149
Unmatched 15843 36
Discarded 0 0

7.1 Balance and estimate

[Link](m_caliper,
stats = "[Link]",
thresholds = c(m = 0.1),
abs = TRUE,
[Link] = "unadjusted",
title = "Balance: caliper matching (0.1 SD)")

12
 Balance: caliper matching (0.1 SD)
distance
re75
re74
age
educ Sample
black Unadjusted

unem74 Adjusted

married
unem75
nodegree
hispan

0 1 2 3
Absolute Mean Differences

matched_cal <- [Link](m_caliper)

model_cal <- lm(re78 ~ treat, data = matched_cal, weights = weights)
coeftest(model_cal, vcov. = vcovHC(model_cal, cluster = ~ subclass))

t test of coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 4708.31 498.89 9.4376 < 0.0000000000000002 ***
treat 1748.42 743.04 2.3531 0.01927 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Is the balance better? Is the estimate closer to $1,794? Try different caliper values (e.g., 0.2,
0.05) and see how the results change.

8 Mahalanobis distance matching

Instead of collapsing all covariates into a single propensity score, we can match directly on
the multivariate distance between covariate vectors. The Mahalanobis distance accounts for
correlations between variables.

13
m_maha <- matchit(treat ~ age + educ + black + hispan + nodegree +
married + re74 + re75 + unem74 + unem75,
data = lalonde,
method = "nearest",
distance = "mahalanobis",
[Link] = "data",
replace = FALSE)

m_maha

A matchit object
- method: 1:1 nearest neighbor matching without replacement
- distance: Mahalanobis
- number of obs.: 16177 (original), 370 (matched)
- target estimand: ATT
- covariates: age, educ, black, hispan, nodegree, married, re74, re75, unem74, unem75

[Link](m_maha,
stats = "[Link]",
thresholds = c(m = 0.1),
abs = TRUE,
[Link] = "unadjusted",
title = "Balance: Mahalanobis distance matching")

Balance: Mahalanobis distance matching

re75

re74

age

educ
Sample
black
Unadjusted
unem74
Adjusted
married

unem75

nodegree

hispan

0 1 2 3
Absolute Mean Differences

14
matched_maha <- [Link](m_maha)
model_maha <- lm(re78 ~ treat, data = matched_maha, weights = weights)
coeftest(model_maha, vcov. = vcovHC(model_maha, cluster = ~ subclass))

t test of coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 4910.69 471.96 10.4050 < 0.0000000000000002 ***
treat 1438.46 747.75 1.9237 0.05516 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

How does this compare to propensity score matching?

9 Full matching

Full matching is an alternative where every unit (treated and control) is placed into a sub-
class, with at least one treated and one control unit per subclass. This avoids discarding any
observations and can improve efficiency. (Note: this method requires the optmatch package.)

m_full <- matchit(treat ~ age + educ + black + hispan + nodegree +

married + re74 + re75 + unem74 + unem75,
data = lalonde,
method = "full",
distance = "glm")

m_full

A matchit object
- method: Optimal full matching
- distance: Propensity score
- estimated with logistic regression
- number of obs.: 16177 (original), 16177 (matched)
- target estimand: ATT
- covariates: age, educ, black, hispan, nodegree, married, re74, re75, unem74, unem75

15
[Link](m_full,
stats = "[Link]",
thresholds = c(m = 0.1),
abs = TRUE,
[Link] = "unadjusted",
title = "Balance: full matching")

Balance: full matching

distance
re75
re74
age
educ Sample
black Unadjusted

unem74 Adjusted

married
unem75
nodegree
hispan

0 1 2 3
Absolute Mean Differences

matched_full <- [Link](m_full)

model_full <- lm(re78 ~ treat, data = matched_full, weights = weights)
coeftest(model_full, vcov. = vcovHC(model_full, cluster = ~ subclass))

t test of coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 5193.60 531.95 9.7632 <0.0000000000000002 ***
treat 1155.54 787.00 1.4683 0.142
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

16
10 Comparing all estimates

Let’s gather our results:

library(broom)

results <- bind_rows(

tidy(model_naive, [Link] = TRUE),
tidy(model_ols, [Link] = TRUE),
tidy(model_nn, [Link] = TRUE),
tidy(model_cal, [Link] = TRUE),
tidy(model_maha, [Link] = TRUE),
tidy(model_full, [Link] = TRUE),
.id = "model"
) |>
filter(term == "treat") |>
mutate(method = c("Naive (no covariates)",
"OLS (with covariates)",
"NN PS matching",
"Caliper matching (0.1 SD)",
"Mahalanobis matching",
"Full matching"))

ggplot(results, aes(x = reorder(method, estimate), y = estimate)) +

annotate("rect",
xmin = -Inf, xmax = Inf,
ymin = 1794 - 100, ymax = 1794 + 100,
fill = "grey80", alpha = 0.5) +
geom_pointrange(aes(ymin = [Link], ymax = [Link])) +
coord_flip() +
labs(x = NULL,
y = "Estimated treatment effect",
title = "Comparing estimation strategies",
caption = "Shaded band = experimental benchmark (~$1,790)") +
theme_minimal()

17
 Comparing estimation strategies

NN PS matching

Caliper matching (0.1 SD)

Mahalanobis matching

Full matching

OLS (with covariates)

Naive (no covariates)

−10000 −5000 0
Estimated treatment effect
Shaded band = experimental benchmark (~$1,790)

18