These notes will examine the incorportion of machine learning methods in classic econometric techniques for estimating causal effects. More specifally, we will focus on estimating treatment effects using matching and instrumental variables. In these estimators (and many others) there is a low-dimensional parameter of interest, such as the average treatment effect, but estimating it requires also estimating a potentially high dimensional nuisance parameter, such as the propensity score. Machine learning methods were developed for prediction with high dimensional data. It is then natural to try to use machine learning for estimating high dimensional nuisance parameters. Care must be taken when doing so though because the flexibility and complexity that make machine learning so good at prediction also pose challenges for inference.

The simplest example of the setting we will analyze is a partially linear model. We have some regressor of interest, $d$, and we want to estimate the effect of $d$ on $y$. We have a rich enough set of controls that we are willing to believe that $\Er[\epsilon|d,x] = 0$. $d_i$ and $y_i$ are scalars, while $x_i$ is a vector. We are not interested in $x$ per se, but we need to include it to avoid omitted variable bias.

Typical applied econometric practice would be to choose some transfrom of $x$, say $X = T(x)$, where $X$ could be some subset of $x$, along with interactions, powers, and so on. Then estimate a linear regression $$y = \theta d + X'\beta + \epsilon$$ and then perhaps also report results for a handful of different choices of $T(x)$.

Some downsides to the typical applied econometric practice include:

• The choice of T is arbitrary, which opens the door to specification searching and p-hacking.

• If $x$ is high dimensional, and $X$ is low dimensional, a poor choice will lead to omitted variable bias. Even if $x$ is low dimensional, if $f(x)$ is poorly approximated by $X’\beta$, there will be omitted variable bias.

In some sense, machine learning can be thought of as a way to choose $T$ is an automated and data-driven way. There will be still be a choice of machine learning method and often tuning parameters for that method, so some arbitrary decisions remain. Hopefully though these decisions have less impact.

You may already be familiar with traditional nonparametric econometric methods like series / sieves and kernels. These share much in common with machine learning. What makes machine learning different that traditional nonparametric methods? Machine learning methods appear to have better predictive performance, and arguably more practical data-driven methods to choose tuning parameters. Machine learning methods can deal with high dimensional $x$, while traditional nonparametric methods focus on situations with low dimensional $x$.

Example : Effect of abortion on crime

Donohue and Levitt (2001) estimate a regression of state crime rates on crime type relevant abortion rates and controls, $$y_{it} = \theta a_{it} + x_{it}'\beta + \delta_i + \gamma_t + \epsilon_{it}.$$ $a_{it}$ is a weighted average of lagged abortion rates in state $i$, with the weight on the $\ell$th lag equal to the fraction of age $\ell$ people who commit a given crime type. The covariates $x$ are the log of lagged prisoners per capita, the log of lagged police per capita, the unemployment rate, per-capita income, the poverty rate, AFDC generosity at time t − 15, a dummy for concealed weapons law, and beer consumption per cap. Alexandre Belloni, Chernozhukov, and Hansen (2014a) reanalyze this setup using lasso to allow a more flexible specification of controls. They allow for many interactions and quadratic terms, leading to 284 controls.

The partially linear and matching models are closely related. If the conditional mean independence assumption of the partially linear model is strengthing to conditional indepence then the partially linear model is a special case of the matching model with constant treatment effects, $y_i(1) - y_i(0) = \theta$. Thus the matching model can be viewed as a generalization of the partially linear model that allows for treatment effect heterogeneity.

All the expectations in these three formulae involve observable data. Thus, we can form an estimate of $\theta$ be replacing the expectations and conditional expectations with appropriate estimators. For example, to use the first formula, we could estimate a logit model for the probability of treatment, $$\hat{\Pr}(d=1|x_i) = \frac{e^{X_i' \hat{\beta}}}{1+e^{X_i'\hat{\beta}}}$$ where, as above, $X$ is a some chosen transformation of $x_i$. Then we simply take an average to estimate $\theta$. $$\hat{\theta} = \frac{1}{n} \sum_{i=1}^n \frac{y_i d_i}{\hat{\Pr}(d=1|x_i)} - \frac{y_i(1-d_i)} {1-\hat{\Pr}(d=1|x_i)}$$ As in the partially linear model, estimating the parameter of interest, $\theta$, requires estimating a potentially high dimensional nuisance parameter, in this case $\hat{\Pr}(d=1|x)$. Similarly, the second expression would require estimating conditional expectations of $y$ as nuisance parameters. The third expression requires estimating both conditional expecations of $y$ and $d$.

The third expression might appear needlessly complicated, but we will see later that it has some desirable properties that will make using it essential when very flexible machine learning estimators for the conditional expectations are used.

The origin of the name “matching” can be seen in the second expression. One way to estimate that expression would be to take each person in the treatment group, find someone with the same (or nearly the same) $x$ in the control group, difference the outcome of this matched pair, and then average over the whole sample. (Actually this gives the average treatment effect on the treated. For the ATE, you would also have to do the same with roles of the groups switched and average all the differences.) When $x$ is multi-dimensional, there is some ambiguity about what it means for two $x$ values to be nearly the same. An important insight of Rosenbaum and Rubin (1983) is that it is sufficient to match on the propensity score, $P(d=1|x)$, instead.

Example: effectiveness of heart catheterization

Connors et al. (1996) use matching to estimate the effectiveness of heart catheterization in critically ill patients. Their dataset contains 5735 patients and 72 covariates. Athey et al. (2017) reanalyze this data using a variety of machine learning methods.

Both the partially linear model and treatment effects model can be extended to situations with endogeneity and instrumental variables.

Most of the remarks about the partially linear model also apply here.

Hartford et al. (2017) estimate a generalization of this model with $y_i = f(d_i, x_i) +\epsilon$ using deep neural networks.

Example : compulsory schooling and earnings

Angrist and Krueger (1991) use quarter of birth as an instrument for years of schooling to estimate the effect of schooling on earnings. Since compulsory schooling laws typically specify a minimum age at which a person can leave school instead of a minimum years of schooling, people born at different times of the year can be required to complete one more or one less year of schooling. Compulsory schooling laws and their effect on attained schooling can vary with state and year. Hence, Angrist and Krueger (1991) considered specifying $g(x,z)$ as all interactions of quarter of birth, state, and year dummies. Having so many instruments leads to statistical problems with 2SLS.

See Abadie (2003).

Belloni et al. (2017) analyze estimation of this model using Lasso and other machine learning methods.

We are following the setup and notation of Chernozhukov et al. (2018). As in the examples, the dimension of $\theta$ is fixed and small. The dimension of $\eta$ is large and might be increasing with sample size. $T$ is some normed vector space.

Example: partially linear model In the partially linear model,

$$y_i = \theta_0 d_i + f_0(x_i) + \epsilon_i$$

we can let $w_i = (y_i, x_i)$ and $\eta = f$. There are a variety of candidates for $\psi$. An obvious (but flawed) one is $\psi(w_i; \theta, \eta) = (y_i - \theta_0 d_i - f_0(x_i))d_i$. With this choice of $\psi$, we have

$$\begin{align*} 0 = & \En[d_i(y_i - \hat{\theta} d_i - \hat{f}(x_i)) ] \\ \hat{\theta} = & \En[d_i^2]^{-1} \En[d_i (y_i - \hat{f}(x_i))] \\ (\hat{\theta} - \theta_0) = & \En[d_i^2]^{-1} \En[d_i \epsilon_i] + \En[d_i^2]^{-1} \En[d_i (f_0(x_i) - \hat{f}(x_i))] \end{align*}$$

The first term of this expression is quite promising. $d_i$ and $\epsilon_i$ are both finite dimensional random variables, so a law of large numbers will apply to $\En[d_i^2]$, and a central limit theorem would apply to $\sqrt{n} \En[d_i \epsilon_i]$. Unfortunately, the second expression is problematic. To accomodate high dimensional $x$ and allow for flexible $f()$, machine learning estimators must introduce some sort of regularization to control variance. This regularization also introduces some bias. The bias generally vanishes, but at a slower than $\sqrt{n}$ rate. Hence

$$\sqrt{n} \En[d_i (f_0(x_i) - \hat{f}(x_i))] \to \infty.$$

To get around this problem, we must modify our estimate of $\theta$. Let $m(x) = \Er[d|x]$ and $\mu(y) = \Er[y|x]$. Let $\hat{m}()$ and $\hat{\mu}()$ be some estimates. Then we can estimate $\theta$ by partialling out:

$$\begin{align*} 0 = & \En[(d_i - \hat{m}(x_i))(y_i - \hat{\mu}(x_i) - \theta (d_i - \hat{m}(x_i)))] \\ \hat{\theta} = & \En[(d_i -\hat{m}(x_i))^2]^{-1} \En[(d_i - \hat{m}(x_i))(y_i - \hat{\mu}(x_i))] \\ (\hat{\theta} - \theta_0) = & \En[(d_i -\hat{m}(x_i))^2]^{-1} \left(\En[(d_i - \hat{m}(x_i))\epsilon_i] + \En[(d_i - \hat{m}(x_i))(\mu(x_i) - \hat{\mu}(x_i))] \right) \\ = & \En[(d_i -\hat{m}(x_i))^2]^{-1} \left( a + b +c + d \right) \end{align*}$$

where

$$\begin{align*} a = & \En[(d_i -m(x_i))\epsilon_i] \\ b = & \En[(m(x_i)-\hat{m}(x_i))\epsilon_i] \\ c = & \En[v_i(\mu(x_i) - \hat{\mu}(x_i))] \\ d = & \En[(m(x_i) - \hat{m}(x_i))(\mu(x_i) - \hat{\mu}(x_i))] \end{align*}$$

with $v_i = d_i - \Er[d_i | x_i]$. The term $a$ is well behaved and $\sqrt{n}a \leadsto N(0,\Sigma)$ under standard conditions. Although terms $b$ and $c$ appear similar to the problematic term in the initial estimator, they are better behaved because $\Er[v|x] = 0$ and $\Er[\epsilon|x] = 0$. This makes it possible, but difficult to show that $\sqrt{n}b \to_p = 0$ and $\sqrt{n} c \to_p = 0$, see e.g. Alexandre Belloni, Chernozhukov, and Hansen (2014a). However, the conditions on $\hat{m}$ and $\hat{\mu}$ needed to show this are slightly restrictive, and appropriate conditions might not be known for all estimators. Chernozhukov et al. (2018) describe a sample splitting modification to $\hat{\theta}$ that allows $\sqrt{n} b$ and $\sqrt{n} c$ to vanish under weaker conditions (essentially the same rate condition as needed for $\sqrt{n} d$ to vanish.)

The last term, $d$, is a considerable improvement upon the first estimator. Instead of involving the error in one estimate, it now involes the product of the error in two estimates. By the Cauchy-Schwarz inequality, $$d \leq \sqrt{\En[(m(x_i) - \hat{m}(x_i))^2]} \sqrt{\En[(\mu(x_i) - \hat{\mu}(x_i))^2]}.$$ So if the estimates of $m$ and $\mu$ converge at rates faster than $n^{-1/4}$, then $\sqrt{n} d \to_p 0$. This $n^{-1/4}$ rate is reached by many machine learning estimators.

Bold references are recommended reading. They are generally shorter and less technical than some of the others. Aspiring econometricians should read much more than just the bold references.

There is considerable overlap among these categories. The papers listed under Neyman orthogonalization all include use of lasso and some include random forests. The papers on lasso all involve some use of orthogonalization.