Example 6: matching using MatchIt

We continue analyzing the lalonde data obtained from the MatchIt package. We estimate the propensity scores using the selected logstic regression model in R Example 5

library("MatchIt")
data("lalonde")

model <- glm(treat ~ re74 + race + married + I(re74^2) + re74:race, data = lalonde, family = "binomial")
eps <- predict(model, type = "response")
lps <- predict(model)

Assess initial imbalance

m.out0 <- matchit(treat ~ age + educ + race + married + 
                   nodegree + re74 + re75, data = lalonde,
                 method = NULL, distance = "glm")
## distance here refer to a method to estimate the propensity score
## No matching is performed at this stage
summary(m.out0)

Call:
matchit(formula = treat ~ age + educ + race + married + nodegree + 
    re74 + re75, data = lalonde, method = NULL, distance = "glm")

Summary of Balance for All Data:
           Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
distance          0.5774        0.1822          1.7941     0.9211    0.3774   0.6444
age              25.8162       28.0303         -0.3094     0.4400    0.0813   0.1577
educ             10.3459       10.2354          0.0550     0.4959    0.0347   0.1114
raceblack         0.8432        0.2028          1.7615          .    0.6404   0.6404
racehispan        0.0595        0.1422         -0.3498          .    0.0827   0.0827
racewhite         0.0973        0.6550         -1.8819          .    0.5577   0.5577
married           0.1892        0.5128         -0.8263          .    0.3236   0.3236
nodegree          0.7081        0.5967          0.2450          .    0.1114   0.1114
re74           2095.5737     5619.2365         -0.7211     0.5181    0.2248   0.4470
re75           1532.0553     2466.4844         -0.2903     0.9563    0.1342   0.2876

Sample Sizes:
          Control Treated
All           429     185
Matched       429     185
Unmatched       0       0
Discarded       0       0

## Love plot before matching
plot(summary(m.out0), abs = F)

## Default thresholds: solid line 0.1, dashed line 0.05

Drop control and treated units that are out of the estimated propenstiy score range

library(ggplot2)
temp.data <- data.frame(lps = lps, treated = as.factor(lalonde$treat))
ggplot(temp.data, aes(x = lps, fill = treated, color = treated)) + 
  geom_histogram(alpha = 0.5, position = "identity") + 
  xlab("Linearized propensity score") 

idx.treated <- which((lalonde$treat == 1) & (lps > max(lps[lalonde$treat == 0])))
idx.control <- which((lalonde$treat == 0) & (lps < min(lps[lalonde$treat == 1])))
idx <- c(idx.treated, idx.control)  ## dropped two treated units and 77 control units
lalonde <- lalonde[-idx, ]
lps <- lps[-idx]

From the histogram, we see a potential issue here, which is that there are not enough controls that can be potentially matched with the treated units for units whose lps > 0.

Matching

Greedy algorithm, covariate balancing is not satisfying

## Greedy algorithm
m.out1 <- matchit(treat ~ age + educ + race + married + 
                   nodegree + re74 + re75, data = lalonde,
                  estimand = "ATT", 
                 method = "nearest", distance = lps)
plot(summary(m.out1), abs = F)

We change to optimal matching, and covariate balancing is still not satisfying

## Optimal matching algorithm
## install additional package
##install.packages("optmatch") and choose the binary version
library("optmatch")

Attaching package: ‘optmatch’

The following object is masked from ‘package:survival’:

    strata

m.out2 <- matchit(treat ~ age + educ + race + married + 
                   nodegree + re74 + re75, data = lalonde,
                  estimand = "ATT", 
                 method = "optimal", distance = lps)
plot(summary(m.out2), abs = F)

You see that even the linearized propensity score are not matched well, which indicates that some treated units are not able to be well matched (as there are not enough controls with lps > 0).

Three strategies to improve balancing:

## Strategy I: require exact match for race, as race seems to be the most unbalanced covariate and is discrete
m.out3 <- matchit(treat ~ age + educ + race + married + 
                   nodegree + re74 + re75, data = lalonde,
                  estimand = "ATT", exact = "race",
                 method = "optimal", distance = lps)

Warning: Fewer control units than treated units in some `exact` strata; not all treated units will get a match.

summary(m.out3)

Call:
matchit(formula = treat ~ age + educ + race + married + nodegree + 
    re74 + re75, data = lalonde, method = "optimal", distance = lps, 
    estimand = "ATT", exact = "race")

Summary of Balance for All Data:
           Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
distance          0.2259       -1.8038          1.9434     0.4496    0.3409   0.6040
age              25.8197       26.4091         -0.0819     0.4847    0.0685   0.1514
educ             10.3279       10.1193          0.1041     0.4802    0.0372   0.1027
raceblack         0.8415        0.2472          1.6276          .    0.5944   0.5944
racehispan        0.0601        0.1733         -0.4762          .    0.1132   0.1132
racewhite         0.0984        0.5795         -1.6158          .    0.4812   0.4812
married           0.1858        0.4205         -0.6033          .    0.2347   0.2347
nodegree          0.7104        0.6335          0.1694          .    0.0769   0.0769
re74           1785.3081     3243.4437         -0.3770     0.7722    0.1575   0.3977
re75           1448.6596     1860.5826         -0.1318     1.3614    0.0915   0.2458

Summary of Balance for Matched Data:
           Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
distance         -0.0701       -0.1097          0.0379     1.0038    0.0095   0.0948
age              26.0086       26.2931         -0.0395     0.4074    0.0929   0.2241
educ             10.4741       10.0431          0.2152     0.3270    0.0517   0.1466
raceblack         0.7500        0.7500          0.0000          .    0.0000   0.0000
racehispan        0.0948        0.0948          0.0000          .    0.0000   0.0000
racewhite         0.1552        0.1552          0.0000          .    0.0000   0.0000
married           0.2328        0.2931         -0.1552          .    0.0603   0.0603
nodegree          0.7069        0.6466          0.1330          .    0.0603   0.0603
re74           2261.7502     2786.3308         -0.1356     0.6877    0.0323   0.1293
re75           2073.3520     1608.2498          0.1488     1.4552    0.0724   0.1552
           Std. Pair Dist.
distance            0.0442
age                 1.1684
educ                1.3341
raceblack           0.0000
racehispan          0.0000
racewhite           0.0000
married             0.1995
nodegree            0.7792
re74                0.5037
re75                0.5826

Sample Sizes:
          Control Treated
All           352     183
Matched       116     116
Unmatched     236      67
Discarded       0       0

plot(summary(m.out3), abs = F)

Strategy II: we can to improve balancing by discard units if the distance between pairs are two large

head(m.out2$match.matrix)

     [,1]     
NSW1 "PSID271"
NSW2 "PSID185"
NSW3 "PSID388"
NSW4 "PSID386"
NSW5 "PSID372"
NSW6 "PSID370"

names(lps) <- rownames(lalonde)
dists <- lps[rownames(m.out2$match.matrix)] - lps[m.out2$match.matrix[,1]]
hist(dists, main = "Histogram of distances between matched pairs")

## We see that some distances are large.

## We set a threshold 1 to make a balance between removing poor matches and retain most of the treated units
ID.treated.notmatched <- names(dists[abs(dists) > 1])
lalonde.new <- lalonde[!(rownames(lalonde) %in% ID.treated.notmatched),]
lps.new <- lps[!(rownames(lalonde) %in% ID.treated.notmatched)]

m.out4 <- matchit(treat ~ age + educ + race + married + 
                   nodegree + re74 + re75, data = lalonde.new,
                  estimand = "ATT", 
                 method = "optimal", distance = lps.new)
plot(summary(m.out4), abs = F)

summary(m.out4)

Call:
matchit(formula = treat ~ age + educ + race + married + nodegree + 
    re74 + re75, data = lalonde.new, method = "optimal", distance = lps.new, 
    estimand = "ATT")

Summary of Balance for All Data:
           Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
distance          0.0178       -1.8038          1.4298     0.6690    0.3048   0.5148
age              25.6727       26.4091         -0.1050     0.4602    0.0689   0.1483
educ             10.2909       10.1193          0.0871     0.4641    0.0379   0.1011
raceblack         0.7364        0.2472          1.1103          .    0.4892   0.4892
racehispan        0.1000        0.1733         -0.2443          .    0.0733   0.0733
racewhite         0.1636        0.5795         -1.1242          .    0.4159   0.4159
married           0.1545        0.4205         -0.7356          .    0.2659   0.2659
nodegree          0.7182        0.6335          0.1882          .    0.0847   0.0847
re74           1209.8659     3243.4437         -0.5896     0.6141    0.2191   0.4818
re75           1187.4235     1860.5826         -0.2158     1.3561    0.1407   0.3301

Summary of Balance for Matched Data:
           Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
distance          0.0178       -0.0923          0.0865     1.0704    0.0261   0.2727
age              25.6727       25.7273         -0.0078     0.4029    0.0914   0.2273
educ             10.2909        9.9455          0.1754     0.3844    0.0478   0.1182
raceblack         0.7364        0.7364          0.0000          .    0.0000   0.0000
racehispan        0.1000        0.0818          0.0606          .    0.0182   0.0182
racewhite         0.1636        0.1818         -0.0491          .    0.0182   0.0182
married           0.1545        0.2455         -0.2515          .    0.0909   0.0909
nodegree          0.7182        0.6545          0.1415          .    0.0636   0.0636
re74           1209.8659     2086.6191         -0.2542     0.6268    0.0901   0.2818
re75           1187.4235     1222.3671         -0.0112     1.3399    0.0224   0.1091
           Std. Pair Dist.
distance            0.0870
age                 1.3358
educ                1.3203
raceblack           0.0000
racehispan          0.0606
racewhite           0.0491
married             0.3018
nodegree            0.9093
re74                0.5885
re75                0.5532

Sample Sizes:
          Control Treated
All           352     110
Matched       110     110
Unmatched     242       0
Discarded       0       0

Covariate balancing improves. However, notice that because we have removed more than 1/3 of of the treated units, the treated population may have changed dramastically. The ATT we estimate here can be biased for the ATT of the original population.

Strategy III: In the software tutorial, it uses another method called “optimal full matching” and reaches extremely good balancing. It assigns every treated and control unit in the sample to one subclass each (Hansen 2004; Stuart and Green 2008). Each subclass contains one treated unit and one or more control units or one control units and one or more treated units. Here it performs better than optimal matching, as full matching allows one control be matched with many treated units.

We can also use matching with replacement here, which does a similar job.

## We match each treated with 2 controls to increase stability (as some controls may be used multiple times)
m.out5 <- matchit(treat ~ age + educ + race + married + 
                   nodegree + re74 + re75, data = lalonde,
                  estimand = "ATT", replace= T, ratio = 2,
                 method = "nearest", distance = lps)
plot(summary(m.out5), abs = F)

Actually, for matching with replacement you can change the propensity score estimation method to the simple logistic regression using all covariates as linear terms, and you will see an even better covariate balancing result.

m.out6 <- matchit(treat ~ age + educ + race + married + 
                   nodegree + re74 + re75, data = lalonde,
                  estimand = "ATT", replace= T, ratio = 2,
                 method = "nearest", distance = "glm") 
plot(summary(m.out6), abs = F)

estimate the ATT

We can continue with either m.out3 or m.out6.

• If we continue with m.out3:

m.data <- match.data(m.out3)
head(m.data)

## Neyman's approach, treating the data as from a pairwise randomized experiment
tau_hat_vec <- sapply(levels(m.data$subclass), function(sc) {
  treated <- m.data$re78[m.data$subclass == sc & m.data$treat == 1]
  control <- m.data$re78[m.data$subclass == sc & m.data$treat == 0]
  return(treated - control)
})

tau_hat <- mean(tau_hat_vec, na.rm = T)
sd_tau_hat <- sd(tau_hat_vec)/sqrt(length(tau_hat_vec))
print(c(tau_hat, sd_tau_hat))

[1] 1461.558 1055.290

  
## 95% CI
print(c(tau_hat- 1.96 * sd_tau_hat, tau_hat + 1.96 * sd_tau_hat))

[1] -606.8106 3529.9267

Neyman’s approach with regression adjustment of the bias

## First run a linear regression on the control (remember to remove the treatment assignment variable!)

fit0 <- lm(re78 ~ age + educ + race + married + nodegree + 
             re74 + re75, data = m.data[m.data$treat == 0, ])

m.data$predicted_y0 <- predict(fit0, m.data)

## Neyman's approach with regression adjustment
tau_hat_vec_adjusted <- sapply(levels(m.data$subclass), function(sc) {
  treated <- m.data$re78[m.data$subclass == sc & m.data$treat == 1]
  control <- m.data$re78[m.data$subclass == sc & m.data$treat == 0]
  
  bias <- m.data$predicted_y0[m.data$subclass == sc & m.data$treat == 1] - 
    m.data$predicted_y0[m.data$subclass == sc & m.data$treat == 0]
  control.adjusted <- control + bias
    
  return(treated - control.adjusted)
})

tau_hat_adjusted <- mean(tau_hat_vec_adjusted)
sd_tau_hat_adjusted <- sd(tau_hat_vec_adjusted)/sqrt(length(tau_hat_vec_adjusted))
print(c(tau_hat_adjusted, sd_tau_hat_adjusted))
  
## 95% CI
print(c(tau_hat_adjusted- 1.96 * sd_tau_hat_adjusted, tau_hat_adjusted + 1.96 * sd_tau_hat_adjusted))

Run outcome regression directly:

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

Results are similar to the regression adjustment of bias.

• If we continue with m.data6

m.data <- match.data(m.out6)
head(m.out6$match.matrix)

## Neyman's approach 
tau_hat_vec <- sapply(rownames(m.out6$match.matrix), function(ID) {
  treated <- m.data[ID,]$re78
  control <- m.data[m.out6$match.matrix[ID,1],]$re78
  return(treated - control)
})


tau_hat <- mean(tau_hat_vec, na.rm = T)

sd_tau_hat <- sd(tau_hat_vec)/sqrt(length(tau_hat_vec))
print(c(tau_hat, sd_tau_hat))
  
## 95% CI
print(c(tau_hat- 1.96 * sd_tau_hat, tau_hat + 1.96 * sd_tau_hat))

Neyman’s approach with regression adjustment of the bias

## First run a linear regression on the control (remember to remove the treatment assignment variable!)

fit0 <- lm(re78 ~ age + educ + race + married + nodegree + 
             re74 + re75, data = m.data[m.data$treat == 0, ])

m.data$predicted_y0 <- predict(fit0, m.data)

## Neyman's approach with regression adjustment
tau_hat_vec_adjusted <- sapply(rownames(m.out6$match.matrix), function(ID) {
  treated <- m.data[ID,]$re78
  control <- m.data[m.out6$match.matrix[ID,1],]$re78
  
  bias <- m.data[ID,]$predicted_y0 - 
    m.data[m.out6$match.matrix[ID,1],]$predicted_y0
  control.adjusted <- control + bias
    
  return(treated - control.adjusted)
})

tau_hat_adjusted <- mean(tau_hat_vec_adjusted)
sd_tau_hat_adjusted <- sd(tau_hat_vec_adjusted)/sqrt(length(tau_hat_vec_adjusted))
print(c(tau_hat_adjusted, sd_tau_hat_adjusted))
  
## 95% CI
print(c(tau_hat_adjusted- 1.96 * sd_tau_hat_adjusted, tau_hat_adjusted + 1.96 * sd_tau_hat_adjusted))

Results are different because:

• m.out3 only estimate the ATT of a subset of treated units
• m.out6 may highly depend on the outcome of a few controls