Chapter 14 - Matching

14.4.1 Distance Matching with Multiple Matching Variables

Copy link

Let’s start with distance measures. How can we take a lot of different matching variables and figure out which two observations are “closest?” We have to boil them down into a single distance measure.

Finally, we can wipe our brow: there’s a single, largely agreed-upon approach to doing this: the Mahalanobis distance.344344 I certainly won’t say it’s the only way to do it. And data scientists are scratching their heads right now because Mahalanobis isn’t anywhere near the most common pick for them, Euclidian (which skips the step of dividing out shared variation - we’ll get there) is. But come on, give a textbook author a victory where they can find one.

The Mahalanobis distance is fairly straightforward. We’ll start with a slightly simplified version of it. First, take each matching variable and divide its value by its standard deviation. Now, each matching variable has a standard deviation of 1. This makes sure that no variable ends up being weighted more heavily just because it’s on a bigger scale.

After we’ve divided by the standard deviation, we can calculate distance. For a given treated observation \(A\) and a given control observation \(B\), the Mahalanobis distance is the sum of the squares of all the differences between \(A\) and \(B\). Then, after you’ve taken the sum, you take the square root. In other words, it’s the sum of squared residuals we’d get if trying to predict the matching variables of \(A\) using the values of \(B\), if the standard deviation of each matching variable was 1. That’s it! That’s the Mahalanobis distance.

Or it’s a simplified version of it, anyway. The one piece I left out is that we’re not just dividing each variable by its standard deviation. Instead, we’re dividing out the whole covariance matrix from the squared values of the variables. If all the matching variables are unrelated to each other, this simplifies to just dividing each variable by its standard deviation. But if they are related to each other, this removes the relationships between the matching variables, making them independent before you match on them.

This requires some matrix algebra to express, but even if you don’t know matrix algebra, all you really need to know is that \(x'x\) can be read as “square all the \(x\)s, then add them up.” Then you can figure out what’s going on here if you squint and remember that taking the inverse (\(^{-1}\)) means “divide by.” The Mahalanobis distance between two observations \(x_1\) and \(x_2\) is:

where \(x_1\) is a vector of all the matching variables for observation 1, similarly for \(x_2\), and \(S\) is the covariance matrix for all the matching variables.

Why is it good for the relationships between different variables to be taken out? Because this keeps us from our matching relying really strongly on some latent characteristic that happens to show up a bunch of times. For example, say you were matching on beard length, level of arm hair, and score on a “masculinity index.” All three of those things are, to one degree or another, related to “being a male.” Without dividing out the covariance, you’d basically be matching on “being a male” three times. Redundant, and the matching process will probably refuse to match any males to any non-males, even if they’re good matches on other characteristics. If we divide out the covariance, we’re still matching on male-ness, but it doesn’t count multiple times.

Returning once again to our credit card debt example, let’s expand to two matching variables: \(BillApril\) and also the \(Age\), in years, of the individuals. Comparing good ol’ treated observation 10,305 and untreated observation 27,281, we already know that the distance between them in their \(BillApril\) is \(|91,630 - 91,609| = 21\). The difference in ages is \(|23 - 37| = -14\).

Let’s calculate Mahalanobis distance first using the simplified approach where we ignore the relationship between the matching variables. The standard deviation of \(BillApril\) is NT$59,554, and so after we divide all the values of \(BillApril\) by the standard deviation, we get a new difference of \(|1.5386 - 1.5382| = .0004\) standard deviations. Not much! Then, the standard deviation of \(Age\) is 9.22, and so we end up with a distance of \(|2.49 - 4.01| = 1.52\) standard deviations.

We have our differences in standard-deviation terms. Square them up and add them together to get \(.0004^2+1.52^2 = 2.310\), and then take the square root - \(\sqrt{2.310} = 1.520\) - to get the Mahalanobis distance of 1.520. We would then compare this Mahalanobis distance to the distance between row 10,305 and all the other untreated rows, ending up in a match with whichever untreated observation leads to the lowest distance.

Let’s repeat that, but properly, taking into account the step where we divide by the covariance matrix. We have a difference of 21 in \(BillApril\) and a difference of -14 in age. We put that next to its covariance matrix and get a Mahalanobis distance of

\[\begin{multline} d(x_1,x_2) = \left(\left[\begin{array}{c} 21 \\ -14 \end{array}\right]' \left[ \begin{array}{cc} 59,544^2 & 26,138 \\ 26,138 & 85 \end{array}\right]^{-1} \left[ \begin{array}{c} 21 \\ -14 \end{array} \right]\right)^{1/2} = \\ \sqrt{.0000257 + 2.312} = 1.521 \end{multline}\]

Okay, not a huge difference from the simplified version in this particular case (1.521 vs. 1.520). But still definitely worth doing. It’s not like we knew it wouldn’t matter before we calculated it both ways.

The Mahalanobis distance has a convenient graphical interpretation, too. We’re squaring stuff, adding it up, and then taking the square root. Ring any bells? Well, for two matching variables that’s just the distance between those two points, as we know from Pythagoras’ theorem. The Mahalanobis distance is the length of the hypotenuse of a triangle with the difference of the first matching variable for the base and the difference of the second matching variable for the height.

Figure 14.6: Mahalanobis Distance with Two Matching Variables

You can see the distance calculation in Figure 14.6. It’s just how long a straight line between the two points is, after you standardize the axes. That’s all! This generalizes to multiple dimensions, too. Add a third matching variable and you’re still drawing a line between two points; it’s just in three dimensions this time. Look around you at any two objects and imagine a line connecting them. That’s the Mahalanobis distance between them if you’re matching on the three spatial dimensions: width, length, and height.

And so on and so forth. Four dimensions is a bit harder to visualize since we don’t see in four dimensions, but it’s the same idea. Add as many dimensions as you like.

This does lead us to one downside of Mahalanobis matching, and distance matching generally - the curse of dimensionality. The more matching variables you add, the “further away” things seem. Think of it this way - imagine a top-down map of an office building. This top-down view has only two dimensions: north/south and east/west. Wanda on the southwest corner of the building is right next to Taylor who’s also on the southwest. A great match! Now move the camera to the side so we can see that same building in three dimensions: north/south, east/west, and up/down. Turns out Wanda’s in the southwest of the 1st floor, but Taylor’s in the southwest of the 22nd floor. No longer a great match.

The curse of dimensionality means that the more matching variables you add, the less likely you are to find an acceptable match for any given treated observation. There are a few ways to deal with this. One approach is to try to limit the number of matching variables, although that means we’re not going to be closing those back doors - no bueno. Another is to just have a whole bunch of observations. Matches being hard to find is okay if you have eight zillion observations to look through to find one. That sounds good, but also, we only have so many observations available.

A third approach is to simply be a little more forgiving about match quality, or at least the value of the caliper/bandwidth, as the number of dimensions goes up. Sometimes people will just try different bandwidths until the amount of matches looks right. One more disciplined approach is to use your normal bandwidth, but divide the Mahalanobis score by the square root of the number of matching variables (⊕Mahony 2014Mahony, Colin. 2014. “Effects of Dimensionality on Distance and Probability Density in Climate Space.” The Seasons Alter.).

Let’s take Broockman to meet Mahalanobis! We’ll be using one-to-one nearest neighbor matching with a Mahalanobis distance measure, calculated using our three matching variables (\(medianhhincom\), \(blackpercent\), \(leg\_democrat\)). We’ll be matching districts with Black legislators to those without (\(leg\_black\)). There are two parts to the code here: the calculation of matching weights, and the estimation of the treatment effect on the outcome variable. I’ll show how to do both. In the case of the Broockman example, the outcome variable is \(responded\), whether the legislator responded to the email. Of course, the actual design was a bit more complex than just looking at the effect of \(leg\_black\) on \(responded\), but we’ll get to that later. For now we’re just showing the relationship between \(leg\_black\) and \(responded\) after matching on \(medianhhincom\), \(blackpercent\), and \(leg\_democrat\).345345 An important note about the code examples in this chapter: I’ll be showing descriptions of how to use these commands, but when it comes to matching, the default options on some important things vary quite a bit from package to package. What standard errors to use? How to handle ties in nearest-neighbor matching, where multiple observations are equally good matches? Should it try to get an average treatment effect or an average treatment on the treated? Don’t expect the different languages to produce the same results. If you want to see what details lead to the different results, and how you could make them consistent, you’ll need to dive into the documentation of these commands.

library(Matching); library(tidyverse)
br <- causaldata::black_politicians

# Outcome
Y <- br %>%
    pull(responded)
# Treatment
D <- br %>%
    pull(leg_black)
# Matching variables
# Note select() is also in the Matching package, so we specify dplyr
X <- br %>%
    dplyr::select(medianhhincom, blackpercent, leg_democrat) %>%
    as.matrix()

# Weight = 2, oddly, denotes Mahalanobis distance
M <- Match(Y, D, X, Weight = 2, caliper = 1)

# See treatment effect estimate
summary(M)

# Get matched data for use elsewhere. Note that this approach will 
# duplicate each observation for each time it was matched
matched_treated <- tibble(id = M$index.treated,
                          weight = M$weights)
matched_control <- tibble(id = M$index.control,
                          weight = M$weights)
matched_sets <- bind_rows(matched_treated,
                          matched_control) 
# Simplify to one row per observation
matched_sets <- matched_sets %>%
                    group_by(id) %>%
                    summarize(weight = sum(weight))
# And bring back to data
matched_br <- br %>%
    mutate(id = row_number()) %>%
    left_join(matched_sets, by = 'id')

# To be used like this! The standard errors here are wrong
lm(responded~leg_black, data = matched_br, weights = weight)

causaldata black_politicians.dta, use clear download

* Create an id variable based on row number,
* which will be used to locate match IDs
g id = _n

* teffects nnmatch does nearest-neighbor matching using first the outcome 
* variable and the matching variables then the treatment variable.
* The generate option creates variables with the match identities
teffects nnmatch (responded medianhhincom blackpercent leg_democrat) ///
    (leg_black), nneighbor(1) metric(mahalanobis) gen(match)
* Note this will produce a treatment effect by itself
* Doing anything else with the match in Stata is pretty tricky
* (or using a caliper). If you want either of those you may want to 
* look into the mahapick or psmatch2 packages
* Although these won't give the most accurate standard errors

import pandas as pd
import numpy as np
# The more-popular matching tools in sklearn
# are more geared towards machine learning than statistical inference
from causalinference.causal import CausalModel
from causaldata import black_politicians
br = black_politicians.load_pandas().data

# Get our outcome, treatment, and matching variables
# We need these as numpy arrays
Y = br['responded'].to_numpy()
D = br['leg_black'].to_numpy()
X = br[['medianhhincom', 'blackpercent', 'leg_democrat']].to_numpy()

# Set up our model
M = CausalModel(Y, D, X)

# Fit, using Mahalanobis distance
M.est_via_matching(weights = 'maha', matches = 1)

print(M.estimates)
# Note it automatically calcultes average treatments on
# average, on treated, and on untreated/control (ATC)

Perhaps Mahalanobis doesn’t go far enough. We can forget about trying to boil down a bunch of distances into one, and really any concerns about one variable being more important than the others, by making the case that all the variables are infinitely important. No mismatch at all. We’re back to exact matching, and specifically coarsened exact matching, which I first mentioned in the section on matching concepts (⊕Iacus, King, and Porro 2012Iacus, Stefano M, Gary King, and Giuseppe Porro. 2012. “Causal Inference Without Balance Checking: Coarsened Exact Matching.” Political Analysis 20 (1): 1–24.).

In coarsened exact matching, something only counts as a match if it exactly matches on each matching variable. The “coarsened” part comes in because, if you have any continuous variables to match on, you need to “coarsen” them first by putting them into bins, rather than matching on exact values. Otherwise you’d never get any matches. What does this coarsening look like? Let’s take \(BillApril\) from the credit card debt example. We might cut the variable up into ten deciles. So for our standby treated observation on row 10,305, instead of matching in that row’s exact value of NT$91,630 (which no other observation shares), we’d recognize that the value of NT$91,630 is between the 80th and 90th percentiles of \(BillApril\), and so just put it in the category of “between NT$63,153 (80th percentile) and NT$112,110 (90th percentile)” (which about 10% of all observations share).346346 That’s a pretty big bin. We’d want to think carefully about whether we think it’s reasonable to say we’ve closed the back door through \(BillApril\) if the way we do it is comparing someone with a bill of NT$64,000 to someone with a bill of NT$112,000. Maybe we need narrower bins for this variable. Construct those bins carefully!

Once you’ve determined how to bin each continuous variable, you look for exact matches. Each treated observation is kept only if it finds at least one exact match, and is dropped otherwise. Each control observation is kept only if it finds at least one exact match as well, and given a weight corresponding to the number of treated-observation matches it has, divided by the number of control observations matched to that treated observation. Then all of that is further multiplied by the total number of matched control observations divided by the total number of matched treatment observations. In other words, \((MyTreatedMatches/MyControlMatches)\) \(\times\) \((TotalControlMatches/TotalTreatedMatches)\). This process ensures, without a doubt, that there are absolutely no differences in the matching variables (after binning) between the treated and control groups.347347 This doesn’t mean there are truly no differences. As mentioned above, any coarsened continuous variable matches exactly on bins, but not on the actual value, so there are differences in the value within bins. There could also be differences in the variables you don’t use for matching. And if the curse of dimensionality encourages you to leave matching variables out when using coarsened exact matching, uh oh!

Coarsened exact matching does require some serious firepower to work, especially once the curse of dimensionality comes into play. Z. Zhao (⊕2004Zhao, Zhong. 2004. “Using Matching to Estimate Treatment Effects: Data Requirements, Matching Metrics, and Monte Carlo Evidence.” Review of Economics and Statistics 86 (1): 91–107.) looks at an example of some fairly standard data, and an attempt to match on twelve variables. This creates about six million cells after interacting all twelve variables. The data set had far fewer than six million observations. Not everyone is going to find a match.

The need for a really big sample size is no joke or something you can skim over, either. Coarsened exact matching, if applied to moderately-sized samples (or any size sample with too many matching variables), can lead to lots of treated observations being dropped. Black, Lalkiya, and Lerner (⊕2020Black, Bernard S, Parth Lalkiya, and Joshua Y Lerner. 2020. “The Trouble with Coarsened Exact Matching.” Northwestern Law & Econ Research Paper Forthcoming.) replicated five coarsened exact matching studies and found that lots of treatment observations ended up getting dropped, which made the treatment effect estimates much noisier, and can also lead the result to be a poor representation of the average treatment effect if certain kinds of treated observations are more likely to find matches than others. One of those studies was Broockman (2013) - uh oh! If you’ve been wondering why the list of matching variables in the Broockman (2013) code examples is so short, it’s because even this list already leads quite a few observations to be dropped with coarsened exact matching, and any more would worsen the problem. Black, Lalkiya, and Lerner (2020) recommend forgetting the method entirely, although I wouldn’t go that far. I would say that you should probably limit its use to huge datasets, and you must check how many treated observations you lose due to not finding any matches. If that’s a lot, you should switch to another method.

The need for a lot of observations, and the relative ease and low computational requirements of calculation,348348 No square roots or logit models to calculate - it matches or it doesn’t. Plus, you don’t need to really match each observation. You can just count the number of treated and untreated observations in each “bin” and create weights on that basis. makes coarsened exact matching fairly popular in the world of big-data data science. There you might hear the exact matches referred to as “doppelgangers.”349349 One other place that coarsened exact matching pops up is in missing-data imputation. Is your data set incomplete, with some data points gone missing? You could look for some exact matches on all your non-missing variables and see what value of the missing variable they have. I discuss missing-data imputation more in Chapter 22. It’s also fairly popular in the world of enormous administrative data sets, where a researcher has managed to get their hands on, say, governmental records on the entire population of a country.350350 Yes, these data sets exist. You often see them for Nordic countries, where you’re really tracked cradle-to-grave. They also exist in more limited forms in other countries. For example there are a few studies that use the entire IRS database of United States tax information. And before you ask, yes of course every researcher is deeply jealous of the people who have access to this data. Would I give a finger for access to this data? Hey, I’ve got nine others.

Let’s look at Klemick, Mason, and Sullivan (⊕2020Klemick, Heather, Henry Mason, and Karen Sullivan. 2020. “Superfund Cleanups and Children’s Lead Exposure.” Journal of Environmental Economics and Management 100: 1022–89.) as one example of coarsened exact matching at work in neither of those contexts. They are interested in the effect of the Environmental Protection Agency’s (EPA) Superfund Cleanup project on the health of children. This was a project where the EPA located areas subject to heavy environmental contamination and then cleaned them up. Did the cleanup process improve childrens’ health?

Their basic research design was to compare the levels of lead in the blood of children in a neighborhood before the cleanup to after, expecting not only that the lead level will drop after cleanup, but that it will improve more in areas very near the site (less than 2 kilometers) than it will in areas a little farther away (2-5 kilometers), which likely weren’t as affected by the pollution when it was there.351351 This design is a preview of what you’ll see in Chapter 18 on difference-in-differences.

As always, there are back doors between distance to a Superfund site and general health indicators, including blood lead levels. The authors close some of these back doors using coarsened exact matching at the neighborhood level, matching neighborhoods on the basis of which Superfund site they’re near, the age of the housing in the area (what percentage of it was built before 1940), the percentage of the population receiving welfare assistance, the percentage who are African American, and the percentage of the population that receives the blood screenings used to measure the blood lead levels used as the dependent variable in the study.

They start with data on about 380,000 people living within 2 kilometers of a Superfund site, and 900,000 living 2-5 kilometers from a Superfund site. Then, the matching leads to a lot of observations being dropped. The sample after matching is about 201,000 from those living within 2 kilometers of a Superfund site, and 353,000 living 2-5 kilometers from a site.

The large drop in sample may be a concern, but we do see considerably reduced differences between those living close to the site or a little farther away. The matching variables look much more similar, of course, but so do other variables that potentially sit on back doors. Percentage Hispanic was not a matching variable, but half of the pre-matching difference disappears after matching. There are similar reductions in difference for traffic density and education levels.

What do they find? Superfund looks like it worked to reduce blood lead levels in children. Depending on specification, they find without doing any matching that lead levels were reduced by 5.1%-5.5%. With matching the results are a bit bigger, from 7.1%-8.3%.

How can we perform coarsened exact matching ourselves? Let’s head back to Broockman (2013), who performed coarsened exact matching in the original paper.

library(cem); library(tidyverse)
br <- causaldata::black_politicians

# Limit to just the relevant variables and omit missings
brcem <- br %>%
    select(responded, leg_black, medianhhincom, 
    blackpercent, leg_democrat) %>%
    na.omit() %>%
    as.data.frame() # Must be a data.frame, not a tibble

# Create breaks. Use quantiles to create quantile cuts or manually for 
# evenly spaced (You can also skip this and let it do it automatically,
# although you MUST do it yourself for binary variables). Be sure
# to include the "edges" (max and min values). So! Six bins each:
inc_bins <- quantile(brcem$medianhhincom, (0:6)/6)

create_even_breaks <- function(x, n) {
    minx <- min(x)
    maxx <- max(x)
    
    return(minx + ((0:n)/n)*(maxx-minx))
}

bp_bins <- create_even_breaks(brcem$blackpercent, 6)

# For binary, we specifically need two even bins
ld_bins <- create_even_breaks(brcem$leg_democrat,2)

# Make a list of bins
allbreaks <- list('medianhhincom' = inc_bins,
                  'blackpercent' = bp_bins,
                  'leg_democrat' = ld_bins)

# Match, being sure not to match on the outcome
# Note the baseline.group is the *treated* group
c <- cem(treatment = 'leg_black', data = brcem,
         baseline.group =  '1',
         drop = 'responded',
         cutpoints = allbreaks,
         keep.all = TRUE)

# Get weights for other purposes
brcem <- brcem %>%
    mutate(cem_weight = c$w)
lm(responded~leg_black, data = brcem, weights = cem_weight)

# Or use their estimation function att. Note there are many options 
# for these functions including logit or machine-learing treatment 
# estimation. Read the docs!
att(c, responded ~ leg_black, data = brcem)

* If necessary: ssc install cem
causaldata black_politicians.dta, use clear download

* Specify our matching variables and treatment. We can also specify 
* how many bins for our variables - (#2) means 2 bins.
* We MUST DO this for binary variables - otherwise
* it will try to split this binary variable into more bins!
cem medianhhincom blackpercent leg_democrat(#2), tr(leg_black)

* Then we can use the cem_weights variable as iweights
reg responded leg_black [iweight = cem_weights]

import pandas as pd
import statsmodels.formula.api as sm
# There is a cem package but it doesn't seem to work that well
# So we will do this by hand
from causaldata import black_politicians
br = black_politicians.load_pandas().data

# Create bins for our continuous matching variables
# cut creates evenly spaced bins
# while qcut cuts based on quantiles
br['inc_bins'] = pd.qcut(br['medianhhincom'], 6)
br['bp_bins'] = pd.qcut(br['blackpercent'], 6)

# Count how many treated and control observations
# are in each bin
treated = br.loc[br['leg_black'] == 1
].groupby(['inc_bins','bp_bins','leg_democrat']
).size().to_frame('treated')
control = br.loc[br['leg_black'] == 0
].groupby(['inc_bins','bp_bins','leg_democrat']
).size().to_frame('control')

# Merge those counts back in
br = br.join(treated, on = ['inc_bins','bp_bins','leg_democrat'])
br = br.join(control, on = ['inc_bins','bp_bins','leg_democrat'])

# For treated obs, weight is 1 if there are any control matches
br['weight'] = 0
br.loc[(br['leg_black'] == 1) & (br['control'] > 0), 'weight'] = 1
# For control obs, weight depends on total number of treated and control
# obs that found matches
totalcontrols = sum(br.loc[br['leg_black']==0]['treated'] > 0)
totaltreated = sum(br.loc[br['leg_black']==1]['control'] > 0)
# Then, control weights are treated/control in the bin,
# times control/treated overall
br['controlweights'] = (br['treated']/br['control']
)*(totalcontrols/totaltreated)
br.loc[(br['leg_black'] == 0), 'weight'] = br['controlweights']

# Now, use the weights to estimate the effect
m = sm.wls(formula = 'responded ~ leg_black',
weights = br['controlweights'],
data = br).fit()

m.summary()

14.4.2 Entropy Balancing

Copy link

There is a newcomer in the world of distance matching that, while it hasn’t been quite as popular in application relative to other options, seems to have some really nice properties. The newcomer is entropy balancing, which grows out of Hainmueller (⊕2012Hainmueller, Jens. 2012. “Entropy Balancing for Causal Effects: A Multivariate Reweighting Method to Produce Balanced Samples in Observational Studies.” Political Analysis 20 (1): 25–46.).

In other matching methods, you aggregate together the matching variables in some way and match on that. Then you hope that the matching you did removed any differences between treatment and control. And as I’ll discuss in the Checking Match Quality section, you look at whether any differences remain and cross your fingers that they’re gone.

The basic concept of entropy balancing is this: instead of aggregating, matching, and hoping, we instead enforce restrictions on the distance between treatment and control. We say “no difference between treatment and control for these moment conditions for these matching variables,” where moment conditions are things like means, variances, skew, and so on. For example “no difference between treatment and control for the mean of \(X_1\).” You can add a lot of those restrictions. Add another restriction that the mean of \(X_2\) must be the same, and \(X_3\), and so on. The descriptive statistics don’t have to be means, either. Entropy balancing lets you require that things like the variance or skewness are the same, too, if you like.352352 The way I think of entropy balancing myself is “sorta like the method of moments but for matching.” That’s a line that will be complete nonsense to most people reading this. But if it’s not nonsense to you, then it took me all of nine words to explain entropy balancing pretty accurately. People do actually use method of moments to generate matching weights in a very similar manner to entropy balancing, too.

Once you have your set of restrictions, entropy balancing gets to work searching for a set of weights that satisfies those restrictions. As long as the restrictions are possible to satisfy, you’ll get them, too.353353 If you can recall your multivariate calculus, it’s really just applying a Lagrange multiplier. That’s just about it! Entropy balancing gives the assured-no-difference result that coarsened exact matching does, but without having to limit the number of matching variables or drop a bunch of treated observations. Why didn’t we do this before?

This section is short for two reasons: first, the intuitive explanation of entropy balancing is quite simple, but to go into any sort of detail would require going a lot more complex. And second, because entropy balancing is not as well-known as other approaches at the time of this writing. Perhaps this is because it requires a fair bit more technical detail than other methods.354354 Also, more of that detail requires trust in your software programmer - what, are you going to code up the weight-searching method? Okay, fine, maybe you are. But for most people it’s a bit more out of their reach than the other methods in this chapter. “I enforced moment conditions” also sounds a bit more like magic than “I found people with similar values of the matching variables” if you’re trying to explain what you’re doing to someone less well-versed in statistics.

But it’s still worth learning about. Perhaps it will become more popular and expected in the future. Perhaps you’ll use it and that will make it more popular and expected in the future. Or at the very least, it can act as a stand-in for the many other matching approaches that didn’t make this book. And there are a few. What, the eight million combinations of options for the matching procedures we have aren’t enough? No! There’s genetic matching, subclassification, optimal matching, etc. etc..

Let’s code up an entropy balance estimation. The examples will use R and Stata. There is a working Python library empirical_calibration that includes entropy balancing, but you’re in for a ride, as even installing it is not as straightforward as other packages. I’ll leave that one out for now, but you can pursue it on your own.355355 See https://github.com/google/empirical_calibration.

library(ebal); library(tidyverse); library(modelsummary)
br <- causaldata::black_politicians

# Outcome
Y <- br %>%
    pull(responded)
# Treatment
D <- br %>%
    pull(leg_black)
# Matching variables
X <- br %>%
    select(medianhhincom, blackpercent, leg_democrat) %>%
    # Add square terms to match variances if we like
    mutate(incsq = medianhhincom^2,
    bpsq = blackpercent^2) %>%
    as.matrix()

eb <- ebalance(D, X)

# Get weights for usage elsewhere
# Noting that this contains only control weights
br_treat <- br %>%
    filter(leg_black == 1) %>%
    mutate(weights = 1)
br_con <- br %>%
    filter(leg_black == 0) %>%
    mutate(weights = eb$w)
br <- bind_rows(br_treat, br_con)

m <- lm(responded ~ leg_black, data = br, weights = weights)
msummary(m, stars = c('*' = .1, '**' = .05, '***' = .01))

* If necessary: ssc install ebalance
causaldata black_politicians.dta, use clear download

* Specify the treatment and matching variables
* And then in targets() specify which moments to match
* 1 for means, 2 for variances, 3 for skew
* Let's do means and variances for our continuous variables
* and just means for our binary matching variable (leg_democrat)
* and store the resulting weights in wt
ebalance leg_black medianhhincom blackpercent leg_democrat, ///
    targets(2 2 1) g(wt)

* Use pweight = wt to adjust estimates
reg responded leg_black [pw = wt]

14.6.1 Conditional Independence Assumption

Copy link

Like with any statistical method, however well-thought out our methods are, they always rely on assumptions. Using matching can relax some assumptions, like linearity (unless we’re using a regression to estimate a propensity score). But that certainly doesn’t make it assumption-free.

Some assumptions carry over from our other work. The main one of those is the conditional independence assumption. Remember in the Regression chapter where I talked about omitted variable bias being the regression lingo for “open back doors”? Well, the conditional independence assumption is matching lingo for “all the back doors are closed.”

In other words, the conditional independence assumption says that the set of matching variables you’ve chosen is enough to close all back doors. The entirety of the relationship between treatment and outcome is either one of the causal front-door paths you want, or is due to a variable that you’ve measured and included as a matching variable.

That’s just new lingo for an assumption we’ve already covered. What else do we need to think about?

14.6.2 Common Support

Copy link

The matching process, no matter how you do it, assumes that there are appropriate control observations to match with. But what if there aren’t? Chaos ensues! This is called a failure of common support.

Why is this a problem? Let’s take an extreme example. You’re interested in the effect of early-release programs on recidivism. These are programs that allow people who have been sentenced to prison time to leave before their sentence is up, often under certain terms and conditions.360360 Parole is an example of an early-release program. You’ve got a group of treated ex-prisoners who were released early, and a group of untreated ex-prisoners who served their full sentence, and information on whether they committed any additional crimes in the ten years after they got out.

You figure that the prisoner’s behavior while in prison might be related both to getting an early release and to their chances of recidivism. So you match on a “behavior score” given by the prison review board that decides whether someone gets early release.

But uh-oh! Turns out that the review board score, which goes from 1 to 10, translates too directly to early release. Everyone with an 8-10 gets early release, and everyone with a 1-7 doesn’t. So when you go looking for matches for your treated observations, you end up trying to find ex-prisoners who got a score of 8-10 but didn’t get early release. And there aren’t any of those. The analysis lacks common support - there simply aren’t any comparable control observations.

Or imagine that instead of trying to select matches, you’re constructing weights. What weight do you give to the 1-7s so that the control group is comparable to the treatment group? Uh, zero I guess. So everything gets a weight of zero and… you again have nothing to compare.

The problem is still there even if it isn’t quite as stark. Maybe all 9-10s get early release while some 8s do and other 8s don’t. We still don’t have any good matches for the 9-10s. Or maybe there are a handful of 9s and 10s that don’t get early release. Now we have matches, but how heavily do we want to rely on that tiny number of matches?

The assumption of common support says that there must be substantial overlap in the distributions of the matching variables comparing the treated and control observations. Or in the context of propensity scores, there must be substantial overlap in the distribution of the propensity score.

Common support is fairly easy to check for, and there are a few ways to do it.

The first approach to checking for common support is simply to look at the distribution of a variable for the treated group against the distribution of that same variable for the untreated group, as in Figure 14.7 (a), which shows the distribution of the percentage Black in district (\(blackpercent\)) matching variable in our Broockman (2013) study.

Figure 14.7: Common Support of the Percentage Black Variable and the Propensity Score in Broockman (2013) Data

Is there good common support for this variable? It’s a bit iffy. It’s not a problem that the distributions are different - we can use matching or re-weighting to scale up the parts of the control distribution until they’re the same as the treated distribution… assuming that those parts of the control distribution exist to be scaled up. If there’s nobody home at that part of the distribution, not much we can do. And there are certainly some regions here where there’s not much to compare to. There aren’t that many districts with a more than 50% Black population but a non-Black legislator. That’s a problem for common support, as many of the Black legislators we do have represent areas with a 50%-80% Black population.

If there are areas with a 50%-80% Black population but a non-Black legislator, then those areas will receive huge weights and increase noise in the estimate. Or if there aren’t any, then those areas with 50%-80% Black population and a Black legislator will have to be dropped. And that’s the bulk of the treated data! The comparison won’t be all that representative of the actual groups with Black legislators.

This same approach of contrasting the distribution works for propensity scores. In Figure 14.7 (b) you can see the distribution of the propensity score for the treated and control groups. As you would expect, the distribution is further to the right for the treated group than the control group.361361 If it weren’t that would be weird, as it would imply that the untreated group was more likely to be treated than the treated group. That’s not a problem. The problem arises if there are big chunks of the distribution for one group (particularly the treated group) where there is no or very very little weight in the control group. Again, we do have some issues here, where there simply aren’t that many districts with non-Black legislators but high propensity scores. In fact, the propensity is so bunched down near 0 for that group that the density in that region is enormous, squashing down the rest of the graph and making it hard to see.

In the context of propensity scores, we may as well take this opportunity to check something else, too - in propensity score estimation we need to assume that the true propensity score is never exactly 0 or 1.362362 Why? It’s easy to see why when we’re weighting - if you’re a treated observation with a 100% probability of treatment, then your weight would be \(1/0\). Oops! Divide by 0. Similarly, an untreated observation with 0% probability of treatment would get \(1/0\). If we see big parts of the distribution right around 0 or 1, this should raise a concern.

Another way to check for common support when you’re selecting matches is simply to see how successful you are at finding matches. If 1% of your treated observations fail to find an acceptable match, that’s pretty good. If 90% of your treated observations fail to find an acceptable match, then you lack common support.

What to do about support? The first step we can take in dealing with support is to avoid trying to match where we don’t have any support. That is, we want to drop observations that don’t match well.

Some matching methods do this for us. If we’re picking matches with a caliper, then bad matches won’t be selected as matches. Observations without common support won’t make it to the matched sample. Done!

Other methods require a little extra work after the fact. Propensity score approaches, especially inverse probability weighting, won’t automatically drop observations outside the range of common support. For this reason it’s common to “trim” the propensity score. For the treated and control groups separately, look for the range of values of the propensity score where the density distribution is zero (or close to zero, below some cutoff value). Then, drop all observations in the other group that have propensity scores in the range that doesn’t exist in your group (⊕Caliendo and Kopeinig 2008Caliendo, Marco, and Sabine Kopeinig. 2008. “Some Practical Guidance for the Implementation of Propensity Score Matching.” Journal of Economic Surveys 22 (1): 31–72.). In other words, limit the data to the range of common support.

Figure 14.8 repeats Figure 14.7 (b), but then trims to the region of common support, requiring at least ten control observations in a given propensity score bin .02 wide.363363 This “minimum to count as common support” approach I have here is arbitrary. I made it up. There’s not really a rule of thumb on how much is enough support if it’s nonzero. In the data we have left, every observation has at least a few observations that have comparable propensity scores. We now know that when we weight to make each group more like the other, we’re not trying to jam in truly incomparable observations. Only the ones we can use to build a comparable group.

Figure 14.8: Propensity Scores After Trimming to Common Support

This is not the only way to trim for common support. A common approach is to trim outliers - any observation with a propensity score in the top or bottom X% gets trimmed. Or instead of working with percentiles, simply trim anything near 0 or 1, maybe anything from \(0-.01\) and \(.99-1\) - remember, we are assuming we don’t have any true 0 or 1 propensities anyway. Another is to just look for any match, rather than a sufficient number of matches in the region. So for example, look at the minimum and maximum propensity score for the treated observations, and get rid of any control observations outside that range.

All of these approaches to enforcing common support entail dropping data. This brings us to the stinging welt of the whole thing - if this process leads to us dropping a lot of observations, especially treated observations, then at best we can say that we’ll be able to generate an estimate of the treatment effect among the individuals that match well, whatever that means, and at worst we can say the matching has failed. Matching relies on the existence of comparable observations, and those might not be around. Maybe a different matching method would work, but more likely matching isn’t going to be able to close your back doors.

14.6.3 Balance

Copy link

The goal of matching is to close back doors through the matching variables. So, uh, did it? Balance is the assumption that our approach to selecting a matched group has closed back doors for the variables we’re interested in. In theory (and in a perfect world), the weights we select should lead to treated and control groups with the exact same values of the back-door variables. No variation between treatment and control for that back-door variable means there’s no back door any more. It has been closed.

“The exact same values” isn’t really going to happen though. So how can we check if we have acceptable levels of balance? And what happens if we don’t?364364 Most of this section also applies to other contexts in which there shouldn’t be any differences between two groups in a set of variables. You could apply these tools to, for example, a randomized experiment, where we’d expect the randomization to lead to balance. Or something like Regression Discontinuity in Chapter 20, where there should be balance for a certain subset of the data. Like anything else that allows for worse matches, bad balance leads to bias. If the balance is bad, the back doors aren’t being closed, so the treatment effect gets contaminated with all those nasty back doors.

Thankfully we don’t have to do much theorizing about whether there’s a problem. We can check in the matched data we have whether there are meaningful differences between the treated and control groups in variables for which there should be no differences - variables that should be balanced.

And if we find that we do have a balance problem, by any of these methods? Then it’s back to the matching procedure. It’s not uncommon for matching-based estimation to be iterative. Match the data, check balance, and if it’s bad, change something about the matching procedure (more variables, different functional form for a propensity score model, a tighter caliper, etc.) and repeat the process.365365 This is one of the selling points for methods like coarsened exact matching and entropy balancing, which by construction avoid balance problems, to the point where the seminal coarsened exact matching paper is titled “Causal Inference Without Balance Checking.” No iterative matching/balance checking!

A common way of checking for balance is a balance table. Balance tables are pretty straightforward. Take a bunch of variables that should have balance. Then, show some of their summary statistics (usually, mean and standard deviation, maybe some other stuff) separately for the treated and control groups. Finally, do a difference-of-means test for each of the variables to see if the treated and control groups differ.366366 If you managed to make your way this far in the book without having taken basic statistics, well, I’m not sure how I feel about that. But in any case, if you did, you might not have done a difference-in-means test. It’s just a test of whether the mean of a variable is different in two groups. You can even do it using regression by regressing the variable on a treatment indicator. Done. You might add some other bells and whistles like a standardized difference in means, where the difference between treated and control groups is divided by their shared standard deviation - fancy! There’s also a variant on standardized difference-in-means called “standardized bias,” a measure you might see on some balance tables, which is standardized difference-in-means with a slightly different standardization.367367 Specifically this is the difference in means divided by a shared variance term, which you get by adding up the variance in the treatment group with the variance in the control group, multiplying the result by .5, and then taking the squared root (see Rosenbaum and Rubin (1985)).

Often, the balance table is created twice: once with the raw, un-matched data to show the extent of the balance problem that matching is to solve, and then again with the matched data to ideally show that the balance problem is gone. What does good balance look like? It looks like no meaningfully large differences in means. And if you’re doing a test, you should find no statistically significant differences in means, or at least no more than you’d expect by chance.368368 If you’ve got a huge long list of variables on your balance table, then it’s likely not a problem if a couple have significant differences, especially if the difference-of-means itself is meaningfully small.

Balance tables are often fairly easy to code. In R, there are many good options in MatchBalance in the Matching package and sumtable in the vtable package, both of which can either be run both without and with matching weights, or with the output from Match, to compare balance before and after matching. In Stata you can follow a teffects analysis with tebalance summarize (although this won’t do significance tests). For a more standard balance table in Stata you can use the balancetable package, adding weights to check the post-matching balance. In Python, if you’ve created a CausalModel() called M using the causalinference package, you can get a balance table using print(M.summary_stats).

We can see an example of a balance table for Broockman (2013), which is a subset of the output produced by MatchBalance in R following the Mahalanobis matching from earlier in this chapter. Importantly, there’s no one standard way to present a balance table, so the format that your software gives you may not match this one. For each variable, before and after matching, it shows the mean of the variable for the treated group, the mean of the variable for the control group, the standardized difference (a.k.a. standardized bias from earlier in this section), and a p-value from a \(t\)-test for the mean being different in the treated and control groups.

Table 14.1: Broockman (2013) Balance Table Before and After Mahalanobis Matching

What do we see here? First, before matching takes place, there’s a serious amount of imbalance. The mean values are both meaningfully and statistically significantly different between the treatment and control groups. How about after matching? Balance looks pretty good. There’s no meaningfully large difference in the means between treatment and control for any of the three variables we’re looking at. For “Legislator is a Democrat” (\(leg\_democrat\) in the code), there’s literally zero difference in the mean.

We do see a statistically significant (at the 95% level) difference in the percentage of the district that is Black after matching, with a \(t\)-statistic p-value below \(.05\). However, this goes along with a meaningfully very close mean value (.515 vs .513). Additionally, to see if we should be concerned, we might want to look at the full distributions and see whether those are different.369369 The full table also contains a Kolmogorov-Smirnov test in addition to a \(t\)-test. This test checks whether two distributions are the same. The p-value for this test for Black percentage is very high.

Figure 14.9: Distribution of a Matching Variable After Matching in the Treated and Control Groups

Balance tables tend to focus on the mean. Is the mean the same in the treated and control groups? Sure, you could use a balance table to check if other things are the same, but you don’t see that a lot.

And that’s a problem! Would you really believe that you’d closed a back door if the treated and control groups had distributions of a matching variable that looked like Figure 14.9? Sure the means are the same. But there’s clearly a difference there.

How can we check for differences between the treated and control groups in the distribution of a variable? Well… I already sorta answered that with Figure 14.9. Just plot the distribution of the variable for both the (matched) treated and control groups, and see how badly they differ. The only difference from Figure 14.9 is that you’d typically overlay the distributions on top of each other so you can more easily see any differences. Overlaying the distributions from Figure 14.9 in Figure 14.10 makes it pretty easy to spot the problem.

Figure 14.10: Overlaid Distribution of a Matching Variable After Matching in the Treated and Control Groups

We wouldn’t necessarily throw out matching for any mismatch between the distributions, but egregious mismatches like the one in Figure 14.10 might be a cause to go back to the drawing board. You can make these graphs for any continuous matching variable, and you can use the tools for discrete distributions from Chapter 3 to do the same for categorical or binary matching variables.

One place that this overlaid density graph is almost mandatory is when you’re working propensity scores. You’ve only got one real matching variable that matters - the propensity score - so you’d better hope the density lines up pretty well.

You can see this in action in Figure 14.11, which shows overlaid density plots for the Broockman study for the matching variable “Black percentage in district,” and for the propensity score, both after applying inverse probability weights (and trimming any propensity scores from \(0-.02\) or from \(.98-1\)).

Figure 14.11: Percentage Black in District and Propensity Score Distributions Before and After Inverse Probability Weighting

How do these look now? They’re… not perfect. With a really good propensity score model, the propensity score distributions for treated and control observations should really be nearly identical. These are a lot better than they were before the matching, but they’re still pretty dissimilar. The same thing goes for percentage Black in district. This might be a good indication to go back to the propensity score model and try adding some things to it to improve the balance we see.

Overlaid density plots are fairly easy to create by hand, but in some cases there are prepackaged versions we can go with.

library(tidyverse)
br <- causaldata::black_politicians

# We can estimate our own propensity score
m <- glm(leg_black ~ medianhhincom + blackpercent + leg_democrat,
    data = br, family = binomial(link = 'logit'))
# Get predicted values
br <- br %>%
    mutate(ps = predict(m, type = 'response'))

# Create IPW weights
br <- br %>%
    mutate(ipw = case_when(
    leg_black == 1 ~ 1/ps,
    leg_black == 0 ~ 1/(1-ps)))

# Density plots for raw data
ggplot(br, aes(x = medianhhincom, color = factor(leg_black))) + 
    geom_density()

# And using our matching weights
ggplot(br, aes(x = medianhhincom, color = factor(leg_black),
    weight = ipw)) + geom_density()

causaldata black_politicians.dta, use clear download

* If we start with teffects...
teffects ipw (responded) (leg_black medianhhincom ///
    blackpercent leg_democrat, logit)

* We can check the balance of a given matching variable
* using tebalance, density which shows both raw and weighted balance
tebalance density medianhhincom

import pandas as pd
import numpy as np
import seaborn as sns
from causaldata import black_politicians
br = black_politicians.load_pandas().data

# Overlaid treatment/control density in raw data
fig1 = sns.kdeplot(data = br, 
x = 'medianhhincom',
hue = 'leg_black',
common_norm = False)
fig1.plot()

# Start the new plot
fig1.get_figure().clf()

# Add weights from any weighting method to check post-matching density
# Here we have the ipw variable from our propensity score matching
# in previous code blocks (make sure you run those first!)
fig2 = sns.kdeplot(data = br, x = 'medianhhincom',
hue = 'leg_black', common_norm = False, weights = 'ipw')
fig2.plot()

The use of a propensity score allows us to test another aspect of balance. In particular, if our sample is properly balanced and the matching has worked to close back doors, then controlling for the propensity score, we should also see no relationship between the matching variables and treatment status, i.e. balance. This is evident from our propensity score causal diagram, reproduced in Figure 14.12. Shut down \(TreatmentPropensity\) and there should be no path from \(A\), \(B\), or \(C\) to \(Treatment\).

Figure 14.12: Causal Diagram Amenable to Propensity Score Matching

Dehejia and Wahba (⊕2002Dehejia, Rajeev H, and Sadek Wahba. 2002. “Propensity Score-Matching Methods for Nonexperimental Causal Studies.” Review of Economics and Statistics 84 (1): 151–61.) offer a way of testing this implication of our diagram using a stratification test. In a stratification test we repeatedly look within limited ranges of the propensity score. By limiting ourselves to a narrow chunk of the propensity score, we are in effect controlling for it. So split the whole range of propensity scores up into chunks and check whether the matching variables and treatment are related within each chunk.370370 How many chunks? They should be narrow enough that the propensity score itself should be balanced within each chunk. If it’s not, go narrower. But not so narrow that there aren’t enough observations within each strata to do a test.

14.7.1 Estimating Mean Differences

Copy link

We’re nearly to the end of the chapter, and we’re just now getting to estimating a treatment effect? Well, yeah. When it comes to matching, most of the work - picking a matching method, doing the matching and making the weights, checking common support and balance, and maybe going back to tweak the matching method if the balance is bad - is already done before we need to estimate anything. But now we’ve done that stuff, and it’s time.371371 Sound the alarm!

With a set of matching weights in hand, how can we estimate the effect of treatment? Get the weighted mean of the outcome for the treated and control groups and compare them. Okay, that’s it, bye.

All right, well, that’s not really it. That is quite a lot of it, though. If all you’re interested in doing is getting the effect of treatment on the mean of the outcome, comparing the weighted means will do it.

One important detail to keep in mind is that while the treatment effect is that simple to estimate, the standard errors on that treatment effect won’t be quite as simple. Calculating the standard error of a weighted mean is easy. But the standard error you calculate won’t remember that the weights had to be estimated and perhaps that some observations were dropped in the matching process. The step of preprocessing the data to match observations and create weights introduces some uncertainty, and so increases the standard errors. Incorporating this uncertainty into your standard errors is important.

There are a few different ways to calculate the appropriate standard error adjustments, most of them having to do with figuring out the variance of the outcome variable conditional on treatment and on the matching variables (or on the weights that come from the matching variables). However, they’re all quite complex, and not really the aim of this book. Until you’re ready to delve into the econometric papers that describe the standard error adjustments you’d want to do,372372 Once again, Caliendo and Kopeinig (2008) pops up as a good place to start. you’re probably going to be using software that does it for you. So the key takeaway here is “unless you’re willing to put in the time to do it right, don’t calculate your own standard errors for treatment effects under matching, let matching-based estimation software do it for you.”

One exception is the use of bootstrap standard errors, as discussed in Chapter 13 on regression. Bootstrap standard errors can be applied to matching estimates. In application of a bootstrap to matching, you first randomly resample the data with replacement (so you get the same sample size you started with, but some observations are replaced with multiple copies of others), then re-perform the matching from scratch, and finally estimate the treatment effect in your bootstrap sample. Repeat this process many, many times, and the standard deviation of the treatment effect across all the bootstrap samples is the standard error of your treatment effect. You can see how this process naturally incorporates any uncertainty we get from the matching process - it allows that uncertainty to feed right into our treatment effect distribution.

Both intuitively (because you are actively incorporating how the matching process affects the sampling distribution) and mechanically (because bootstrap standard errors are fairly simple to do yourself by hand, it takes only a little coding if it’s not already in the command you’re using), bootstrap standard errors can be a great option for matching estimates. That said, bootstrap standard errors can only be used with the “constructing a weighted matched sample” approach to matching. They simply don’t work properly for the “selecting matches” approach, because the sharp in/out decisions that process makes don’t translate into the bootstrap having the proper sampling distribution (⊕Abadie and Imbens 2008Abadie, Alberto, and Guido W Imbens. 2008. “On the Failure of the Bootstrap for Matching Estimators.” Econometrica 76 (6): 1537–57.). So if you’re using a weighting approach, go for it. If not, either go with an appropriate standard error adjustment in your software, or if you’re feeling real spicy you can go find one of the bootstrap variants designed to work properly with selecting-matches approaches.

Bootstrapping a matching process can be a little tricky to code up if you’re not using a command that offers it already, since you aren’t just bootstrapping a single command (or asking for the “bootstrap standard error” option in a regression). You need to bootstrap the whole matching process from start to finish. I offer some coding examples later in the chapter in application to doubly robust estimation, and also in Chapter 15. Slight tweaks on this code should help you provide bootstrap standard errors for all kinds of matching processes.

14.7.2 Combining Matching and Regression

Copy link

In many cases we want to do something a little more complex than comparing means between the treated and control groups. For example, maybe our research design depends on a specific functional form of the relationship between treatment and outcome, which would be difficult to capture with matching but easy with regression. But we still want to close back doors by matching.373373 This will come back in Chapter 18 on difference-in-differences, which relies on comparing a treated and control group (which might be matched), and pinpointing a particular discontinuity in the relationship between outcome and time. What then? We can combine our matching weights with a regression model to do regression adjustment.

The approach here is pretty straightforward. We already have a set of selected matched observations, or a set of matching weights, from our matching procedure. We can apply that subsample, or those weights, to our regression model, just as we did with sample weights in Chapter 13.

In fact, we’ve already done this procedure - the basic difference in weighted means we talked about before can be achieved by applying matching weights to a regression with only treatment as a predictor, i.e. \(Outcome = \beta_0 + \beta_1Treatment + \varepsilon\). The estimate of the effect would be \(\hat{\beta}_1\). This is the method I’ve been using in the code examples. It’s not a big jump to start adding other things to that regression.

Using weights isn’t the only way to go (and may not be ideal if you already have other plans for those weights, like survey sample weights). Another approach to using matching with regression is to think carefully about what the propensity score is - it’s a variable that, if controlled for, should block all back doors from \(Treatment\) to \(Outcome\). So… control for it! The propensity score (or the inverse probability weight based on the propensity score) can be added to regression as a control.

In fact, now that we’re at this section we can finally come back around and run the analysis that Broockman (2013) actually ran. He didn’t just do matching and estimate a treatment effect from that. No! He took those matching weights and used them in a regression that let him not just see the difference in response rates between Black legislators and others, but see whether the difference between in-district vs. out-of-district response rates for these letters from Black people was different between Black legislators and others. That calls for an interaction term in a regression. Not to mention we can add plenty of other controls.

library(cem); library(tidyverse); library(modelsummary)
br <- causaldata::black_politicians

# This copies the CEM code from the CEM section
# See that section's code for comments and notes

# Limit to just the relevant variables and omit missings
# (of which there are none in this data)
brcem <- br %>%
    select(responded, leg_black, medianhhincom, 
    blackpercent, leg_democrat) %>%
    na.omit() %>%
    as.data.frame() 

# Create breaks
inc_bins <- quantile(brcem$medianhhincom, (0:6)/6)

create_even_breaks <- function(x, n) {
    minx <- min(x)
    maxx <- max(x)
    return(minx + ((0:n)/n)*(maxx-minx))
}

bp_bins <- create_even_breaks(brcem$blackpercent, 6)
ld_bins <- create_even_breaks(brcem$leg_democrat,2)

allbreaks <- list('medianhhincom' = inc_bins,
    'blackpercent' = bp_bins,
    'leg_democrat' = ld_bins)

c <- cem(treatment = 'leg_black', data = brcem,
    baseline.group =  '1', drop = 'responded',
    cutpoints = allbreaks, keep.all = TRUE)

# Get weights for other purposes.  Note this exact code only 
# works because we didn't have to drop any NAs. If we did, 
# lining things up would be trickier
br <- br %>%
    mutate(cem_weight = c$w)
m1 <- lm(responded~leg_black*treat_out + nonblacknonwhite +
    black_medianhh + white_medianhh + statessquireindex + 
    totalpop + urbanpercent, data = br, weights = cem_weight)
msummary(m1, stars = c('*' = .1, '**' = .05, '***' = .01))

causaldata black_politicians.dta, use clear download

* Use any matching method to get weights
* We'll use coarsened exact matching like in the original
cem medianhhincom blackpercent leg_democrat(#2), tr(leg_black)

reg responded i.treat_out##i.leg_black nonblacknonwhite leg_democrat ///
    leg_senator south blackpercent black_medianhh white_medianhh ///
    statessquire totalpop urbanpercent [iweight=cem_weights]

import pandas as pd
import statsmodels.formula.api as sm
from causaldata import black_politicians
br = black_politicians.load_pandas().data

# This copies the CEM code from the CEM section
# See that section's code for comments and notes
br['inc_bins'] = pd.qcut(br['medianhhincom'], 6)
br['bp_bins'] = pd.qcut(br['blackpercent'], 6)
treated = br.loc[br['leg_black'] == 1
].groupby(['inc_bins','bp_bins','leg_democrat']
).size().to_frame('treated')
control = br.loc[br['leg_black'] == 0
].groupby(['inc_bins','bp_bins','leg_democrat']
).size().to_frame('control')

# Merge the counts back in
br = br.join(treated, on = ['inc_bins','bp_bins','leg_democrat'])
br = br.join(control, on = ['inc_bins','bp_bins','leg_democrat'])

# Create weights
br['weight'] = 0
br.loc[(br['leg_black'] == 1) & (br['control'] > 0), 'weight'] = 1
totalcontrols = sum(br.loc[br['leg_black']==0]['treated'] > 0)
totaltreated = sum(br.loc[br['leg_black']==1]['control'] > 0)
br['controlweights'] = (br['treated']/br['control']
)*(totalcontrols/totaltreated)
br.loc[(br['leg_black'] == 0), 'weight'] = br['controlweights']

# Now, use the weights to estimate the effect
m = sm.wls(formula = '''responded ~ leg_black*treat_out + 
nonblacknonwhite + black_medianhh + white_medianhh + 
statessquireindex + totalpop + urbanpercent''',
weights = br['weight'],
data = br).fit()
m.summary()

From these regressions we get a coefficient of -.142 on \(leg\_black\), -.351 on \(treat\_out\), and .229 on \(leg_black\times treat\_out\). Recalling what we learned about interpreting a model with interaction terms from Chapter 13, we know that this means that Black legislators were 14.2 percentage points less likely than White legislators to respond to in-district emails, out-of-district emails were 35.1 percentage points less likely to get a response than in-district emails from White legislators, but the in-district vs. out-of-district gap is 22.9 percentage points smaller for Black legislators than White ones, which is evidence in favor of the idea that Black legislators have an intrinsic motivation to help Black people.

If we play our cards right, regression and matching can combine to be more than the sum of their parts. Doubly robust estimation is a way of combining regression and matching that works even if there’s something wrong with the matching or something wrong with the regression - as long as it’s not both.

What do I mean by “something wrong?” It doesn’t cover everything - you still need to, you know, actually measure for and include the variables that allow you to close back doors.374374 Failing to do this would make neither method work - remember, at least one needs to work. Doubly robust estimation won’t let you get around that. But regression and matching both rely on assumptions. For example, perhaps your regression model is misspecified - you’ve got all the right variables in there, but maybe the functional form isn’t quite right. The same thing can happen when specifying the model you use to estimate your propensity score.

Doubly robust estimation is any method that still works as long as one of those two models is properly specified. The other one can be wrong.

Double-robustness is a thing that some estimators have rather than being a specific estimator. So there are actually lots of methods that are doubly robust.375375 How do we know if a particular estimator is doubly robust? Someone writes a paper with a lot of equations proving that it is. There are a lot of these papers. I could list a dozen right here without even thinking too hard. But we don’t have all day. So we’ll look at three, the first one in detail.

The first estimator we’ll look at is the one described in Wooldridge (2010, Chapter 21.3.4). This one is appealing because the process of doing it is almost like you’d guess doubly robust estimation works on your own if you’d just heard about it and had to come up with it yourself. And it actually does work!

The steps are fairly straightforward:

1. Estimate the propensity score \(p\) for each observation using your matching variables
2. Use the propensity score to produce the inverse probability weights: \(1/p\) for treated observations and \(1/(1-p)\) for untreated observations
3. Using only treated observations, estimate your regression model, using the matching variables as predictors (and perhaps some other predictors too, depending on what you’re doing), and inverse probability weights
4. Repeat step (3) but using only untreated observations
5. Use the models from steps (3) and (4) to produce “treated” and “untreated” predicted observations for the whole sample
6. Compare the “treated” means to the “untreated” means to get the estimate of the causal effect
7. To get standard errors, use bootstrap standard errors, or the heteroskedasticity-robust errors described in Wooldridge (2010) (not just the same as asking for regular heteroskedasticity-robust standard errors)

While there are a lot of steps here, none of them are particularly complex, and all are fairly easy to do by hand. In fact, let’s work through the code for this approach.

library(boot); library(tidyverse)
br <- causaldata::black_politicians

# Function to do IPW estimation with regression adjustment
ipwra <- function(br, index = 1:nrow(br)) {
    # Apply bootstrap index
    br <- br %>% slice(index)
    
    # estimate and predict propensity score
    m <- glm(leg_black ~ medianhhincom + blackpercent + leg_democrat,
             data = br, family = binomial(link = 'logit'))
    br <- br %>%
        mutate(ps = predict(m, type = 'response'))
    
    # Trim control observations outside of treated PS range
    minps <- br %>%
        filter(leg_black == 1) %>%
        pull(ps) %>%
        min(na.rm = TRUE)
    maxps <- br %>%
        filter(leg_black == 1) %>%
        pull(ps) %>%
        max(na.rm = TRUE)
    br <- br %>%
        filter(ps >= minps & ps <= maxps)
    
    # Create IPW weights
    br <- br %>%
        mutate(ipw = case_when(
        leg_black == 1 ~ 1/ps,
        leg_black == 0 ~ 1/(1-ps)))
    
    # Estimate difference
    w_means <- br %>% 
        group_by(leg_black) %>%
        summarize(m = weighted.mean(responded, w = ipw)) %>%
        arrange(leg_black)
    
    return(w_means$m[2] - w_means$m[1])
}


b <- boot(br, ipwra, R = 200)
# See estimate and standard error
b

causaldata black_politicians.dta, use clear download

* Wrap this in a program so we can bootstrap
capture program drop ipwra
program def ipwra, rclass
quietly{    
    * Estimate propensity
    logit leg_black medianhhincom blackpercent leg_democrat
    * Get predicted values for our propensity score
    predict ps, pr
    
    * "Trim" control observations outside of 
    * treated propensity score range
    summ ps if leg_black
    replace ps = . if ps < r(min) | ps > r(max)
    
    * Get inverse probability weights
    g ipw = 1/ps if leg_black
    replace ipw = 1/(1-ps) if !leg_black
    
    * Regress treated and nontreated separately, 
    * then predict for whole sample
    regress responded medianhhincom blackpercent ///
        leg_democrat [pw = ipw] if !leg_black
    predict untreated_outcome
    regress responded medianhhincom blackpercent ///
        leg_democrat [pw = ipw] if leg_black
    predict treated_outcome
    
    * Compare means
    summarize untreated_outcome
    local unt_mean = r(mean)
    summarize treated_outcome
    local tr_mean = r(mean)
    
    local diff = `tr_mean' - `unt_mean'
}

drop ps ipw untreated_outcome treated_outcome
return scalar diff = `diff'
end

* And a version for bootstrapping
capture program drop ipwra_b
program def ipwra_b, rclass
    preserve
    * Create bootstrap sample
    bsample
    * Estimate
    ipwra
    * Bring back original data
    restore
    return scalar diff = r(diff)
end

* Run once with original data
ipwra
local effect = r(diff)
* Then boostrap
simulate diff = r(diff), reps(2000): ipwra_b    
local N = _N
bstat, stat(`effect') n(`N')

* Or skip all that and just use teffects ipwra (no trim here though)
causaldata black_politicians.dta, use clear download
teffects ipwra (responded medianhhincom blackpercent leg_democrat) ///
    (leg_black medianhhincom blackpercent leg_democrat, logit)

import pandas as pd
import numpy as np
import statsmodels.formula.api as sm
from causaldata import black_politicians
br = black_politicians.load_pandas().data

# As mentioned, est_via_weighting from causalinfernece already does 
# doubly robust estimation but it's a different kind! 
# Let's do Wooldridge here.
def ipwra(br):
    # Estimate propensity
    m = sm.logit('leg_black ~ medianhhincom + blackpercent + leg_democrat', 
    data = br)
    # Get fitted values and turn them into probabilities
    m = m.fit().fittedvalues
    br['ps'] = np.exp(m)/(1+np.exp(m))
    
    # Trim observations outside of treated range
    minrange = np.min(br.loc[br['leg_black']==1, 'ps'])
    maxrange = np.max(br.loc[br['leg_black']==1, 'ps'])
    br = br.loc[(br['ps'] >= minrange) & (br['ps'] <= maxrange)]
    
    # Get inverse probability score
    br['ipw'] = br['leg_black']*(1/br['ps']
    ) + (1-br['leg_black'])*(1/(1-br['ps']))
    
    # Regress treated and nontreated separately, 
    # then predict for whole sample 
    mtreat = sm.wls('''responded ~ medianhhincom + 
    blackpercent + leg_democrat''',
    weights =  br.loc[br['leg_black'] == 1, 'ipw'], 
    data = br.loc[br['leg_black'] == 1])
    mcontrol = sm.ols('''responded ~ medianhhincom + 
    blackpercent + leg_democrat''',
    weights =  br.loc[br['leg_black'] == 0, 'ipw'],
    data = br.loc[br['leg_black'] == 0])
    
    treat_predict = mtreat.fit().predict(exog = br)
    con_predict = mcontrol.fit().predict(exog = br)
    
    # Compare means
    diff = np.mean(treat_predict) - np.mean(con_predict)
    return diff

# And a wrapper function to bootstrap
def ipwra_b(br):
    n = br.shape[0]
    br = br.iloc[np.random.randint(n, size=n)]
    diff = ipwra(br)
    return diff

# Run once on the original data to get our estimate
est = ipwra(br)

# And then a bunch of times to get the sampling distribution
dist = [ipwra_b(br) for i in range(0,2000)]

# Our estimate
est
# and its standard error
np.std(dist)

How about the other two methods I mentioned? I’ll talk about these briefly. Just a tour of the space.

One is the Augmented Inverse Probability Weighted Estimator (AIPWE) (⊕Tan 2010Tan, Zhiqiang. 2010. “Bounded, Efficient and Doubly Robust Estimation with Inverse Weighting.” Biometrika 97 (3): 661–82.). Like the Wooldridge approach, AIPWE has both a matching model to estimate, with treatment as a dependent variable, and a regression model, with the outcome as a dependent variable. The basic idea of AIPWE is that it includes an augmentation term (surprise) in the regression model that adjusts the results based on the degree of misspecification in the matching model.

AIPWE has some nice properties. It does okay in some settings even with mild misspecification in both models (although again, you would prefer at least one being correct). Plus, if it turns out that there isn’t any misspecification in the model for treatment, AIPWE will produce lower standard errors in the outcome model. Neat!

The third and final method can be discussed briefly because I’ve already discussed it. Entropy Balancing, which showed up back in the Distance Matching section, is itself doubly robust (⊕Q. Zhao and Percival 2017Zhao, Qingyuan, and Daniel Percival. 2017. “Entropy Balancing Is Doubly Robust.” Journal of Causal Inference 5 (1).). Despite not directly including a regression, entropy balancing’s focus on removing differences in means between treatment and control turns out to work a lot like regression’s focus on removing differences in means. Double robustness!

14.8.1 Getting Other Treatment Effects from Matching

Copy link

Ah, but I have buried the lede.380380 Yes, that’s how it’s spelled. The average treatment on the treated is not the only thing that matching can give you. We can get the average treatment effect, or even the average treatment on the untreated.381381 Granted, subject to the distribution-weighting thing I just talked about. How can we do that?

It’s actually fairly straightforward. If you want average treatment on the treated, as we’ve been getting, then construct a control group that is as similar to the treated group as possible.

Want average treatment on the untreated? Simple! Just reverse things. Instead of matching control observations to treated observations, instead match treated observations to control observations. Same procedure, but reverse which group is which. The result gives you a treated group that is as much like the untreated group as possible. Comparing that to the untreated group compares what you actually got from no treatment against what we think that group would have gotten if treated. Average treatment on the untreated!382382 This does all assume that the variables that predict variation in the treatment effect are included as matching variables. And if you’re using a propensity score to match, it further assumes that the variables that predict variation in the treatment effect are the same ones that predict treatment itself.

From here you can get to the average treatment effect fairly easily. If the average among your treated group is the average treatment effect on the treated (ATT), and the average among your untreated group is the average treatment effect on the untreated (ATUT), then what’s the overall average, also known as the average treatment effect? Well, if a proportion \(p\) of people are getting treated, then it’s just \(ATE = (p\times ATT) + ((1-p)\times ATUT)\). You’re just averaging the average from each “side,” weighted by how big each side is, to get the overall average.

The use of matching to get average treatment effects is fairly common, and sometimes is even the default option in software commands designed to do treatment effects by matching. Stata’s teffects is an example of this - the teffects examples in this chapter produce average treatment effects by default (and all of them have an atet option for getting average treatment on the treated).

How does this apply to inverse probability weighting? After all, when you’re doing inverse probability weighting you aren’t really matching to anything, you’re just weighting on the basis of the inverse propensity weight.

It turns out that the specific way you choose to turn the propensity score into an inverse propensity weight determines the kind of treatment effect average you get. Weighting to make the control group look like the treatment? Average treatment on the treated! Doing the reverse? Average treatment on the untreated! Doing both? Average treatment effect!

In fact, the approach to inverse probability weighting I’ve been talking about this whole time - with a weight of \(1/p\) for the treated group and \(1/(1-p)\) for the untreated group, where \(p\) is the propensity score - serves to both upweight the members of the untreated group that are most like the treated group and upweight the members of the treated group that are most like the untreated group. This actually gives us the average treatment effect, not the average treatment on the treated.

What if we do want one of those other averages? We just use different weights.

What if we want average treatment on the treated? Intuitively here we want to make the untreated group like the treated group. So first we give all of the treated group the same weight of 1. We don’t need to adjust them to be like anyone else - we want to adjust everyone else to be like them so the others can serve as a comparison. Sounds like average treatment on the treated to me.

The treated group gets a weight of 1. How about the untreated group? Instead of a weight of \(1/(1-p)\), they get \(p/(1-p)\). The \(1/(1-p)\), which we did before, sort of brings the representation of the untreated group back to neutral, where everyone is counted equally. But we don’t want equal. We want the untreated group to be like the treated group, which means weighting them even more based on how like the treated group they are. And thus \(p/(1-p)\), which is more heavily responsive to \(p\) than \(1/(1-p)\) is.

From there, average treatment on the untreated is pretty easy to figure out, and works for largely the same reason. The untreated group in this case gets a weight of 1, and the treated group, instead of the \(1/p\) that gave us the ATE, gets \((1-p)/p\), which makes things match the untreated group properly and gives us the ATUT we want.