Double/Debiased/Neyman Machine Learning of Treatment Effects

TL;DR

This paper discusses a double/debiased machine learning approach for valid inference on treatment effects, utilizing Neyman-orthogonal scores and cross-fitting, especially in high-dimensional observational data.

Contribution

It applies the double/debiased machine learning method to estimate average treatment effects and effects on the treated, demonstrating its practical use in observational studies.

Findings

Effective estimation of ATE and ATTE using high-dimensional data.

Validation of the double/debiased ML approach for treatment effect inference.

Enhanced accuracy in treatment effect estimation with modern machine learning methods.

Abstract

Chernozhukov, Chetverikov, Demirer, Duflo, Hansen, and Newey (2016) provide a generic double/de-biased machine learning (DML) approach for obtaining valid inferential statements about focal parameters, using Neyman-orthogonal scores and cross-fitting, in settings where nuisance parameters are estimated using a new generation of nonparametric fitting methods for high-dimensional data, called machine learning methods. In this note, we illustrate the application of this method in the context of estimating average treatment effects (ATE) and average treatment effects on the treated (ATTE) using observational data. A more general discussion and references to the existing literature are available in Chernozhukov, Chetverikov, Demirer, Duflo, Hansen, and Newey (2016).

Tables

Table 1. Table 1. Estimated Mean ATE of 401(k) Eligibility on Net Financial Assets

	Lasso	Reg. Tree	Random Forest	Boosting	Neural Net.	Ensemble	Best
A. Interactive Model
ATE (2 fold)	6331	7581	7966	7826	7805	7617	7800
	(2712)	(1374)	(1549)	(1345)	(1688)	(1299)	(1325)
ATE (5 fold)	6964	8023	8104	7699	7772	7658	7890
	(1654)	(1311)	(1364)	(1223)	(1324)	(1204)	(1198)
B. Partially Linear Model
ATE (2 fold)	7718	8745	9180	8768	9040	9043	9106
	(1796)	(1488)	(1526)	(1451)	(1494)	(1432)	(1430)
ATE (5 fold)	8182	8913	9248	9092	9038	9186	9214
	(1578)	(1440)	(1402)	(1380)	(1394)	(1381)	(1361)
Notes: Estimated Mean ATE and standard errors (in parentheses) from a linear model (Panel B) and heterogeneous effect model (Panel A) based on orthogonal estimating equations. Column labels denote the method used to estimate nuisance functions. Further details about the methods are provided in the main text.

Table 2. Table 2. Estimated Median ATE of 401(k) Eligibility on Net Financial Assets

	Lasso	Reg. Tree	Random Forest	Boosting	Neural Net.	Ensemble	Best
A. Interactive Model
ATE (2 fold)	6725	7557	8034	7820	7800	7620	7800
	(1612)	(1283)	(1400)	(1199)	(1474)	(1198)	(1185)
ATE (5 fold)	7133	8046	8099	7690	7795	7668	7876
	(1420)	(1242)	(1296)	(1179)	(1290)	(1180)	(1149)
B. Partially Linear Model
ATE (2 fold)	7707	8770	9204	8746	9104	9061	9129
	(1785)	(1424)	(1392)	(1391)	(1388)	(1343)	(1342)
ATE (5 fold)	8202	8894	9252	9089	9065	9199	9232
	(1581)	(1440)	(1400)	(1378)	(1393)	(1379)	(1359)
Notes: Estimated Median ATE and standard errors (in parentheses) from a linear model (Panel B) and heterogeneous effect model (Panel A) based on orthogonal estimating equations. Column labels denote the method used to estimate nuisance functions. Further details about the methods are provided in the main text.

Table 3. Table 3. Estimated Mean ATE of Cash Bonus on Unemployment Duration

	Lasso	Reg. Tree	Random Forest	Boosting	Neural Net.	Ensemble	Best
A. Interactive Model
ATE (2 fold)	-0.081	-0.084	-0.072	-0.078	-0.073	-0.079	-0.078
	(0.036)	(0.037)	(0.042)	(0.036)	(0.041)	(0.036)	(0.036)
ATE (5 fold)	-0.081	-0.084	-0.070	-0.076	-0.072	-0.079	-0.076
	(0.036)	(0.037)	(0.040)	(0.036)	(0.038)	(0.036)	(0.036)
B. Partially Linear Model
ATE (2 fold)	-0.081	-0.083	-0.076	-0.076	-0.073	-0.076	-0.076
	(0.036)	(0.037)	(0.037)	(0.036)	(0.036)	(0.036)	(0.036)
ATE (5 fold)	-0.080	-0.084	-0.075	-0.075	-0.074	-0.075	-0.075
	(0.036)	(0.037)	(0.036)	(0.036)	(0.036)	(0.036)	(0.036)
Notes: Estimated Mean ATE and standard errors (in parentheses) from a linear model (Panel B) and heterogeneous effect model (Panel A) based on orthogonal estimating equations. Column labels denote the method used to estimate nuisance functions. Further details about the methods are provided in the main text.

Table 4. Table 4. Estimated Median ATE of Cash Bonus on Unemployment Duration

	Lasso	Reg. Tree	Random Forest	Boosting	Neural Net.	Ensemble	Best
A. Interactive Model
ATE (2 fold)	-0.081	-0.084	-0.073	-0.078	-0.074	-0.079	-0.078
	(0.036)	(0.036)	(0.041)	(0.036)	(0.039)	(0.036)	(0.036)
ATE (5 fold)	-0.081	-0.085	-0.069	-0.076	-0.072	-0.079	-0.076
	(0.036)	(0.037)	(0.039)	(0.036)	(0.038)	(0.036)	(0.036)
B. Partially Linear Model
ATE (2 fold)	-0.081	-0.084	-0.077	-0.076	-0.074	-0.076	-0.076
	(0.036)	(0.036)	(0.036)	(0.036)	(0.036)	(0.036)	(0.036)
ATE (5 fold)	-0.079	-0.084	-0.076	-0.075	-0.073	-0.075	-0.075
	(0.036)	(0.037)	(0.036)	(0.035)	(0.036)	(0.035)	(0.035)
Notes: Estimated Median ATE and standard errors (in parentheses) from a linear model (Panel B) and heterogeneous effect model (Panel A) based on orthogonal estimating equations. Column labels denote the method used to estimate nuisance functions. Further details about the methods are provided in the main text.

Equations

$$Y=g_0(D,Z)+\zeta ,$$

$$D=m_0(Z)+\nu ,$$

$${\theta }_0=\text{E}[g_0(1,Z)-g_0(0,Z)],$$

$${\theta }_0=\text{E}[g_0(1,Z)-g_0(0,Z)|D=1].$$

$$\text{E}\psi (W;{\theta }_0,{\eta }_0)=0,$$

$$\displaystyle\partial_{\eta}\text{E}\psi(W;\theta_{0},\eta)\Big{|}_{\eta=\eta_{0}}=0,$$

$$\begin{array}{rl}\psi (W;\theta ,\eta )& :=(g(1,Z)-g(0,Z))+\frac{D(Y-g(1,Z))}{m(Z)}-\frac{(1-D)(Y-g(0,Z))}{1-m(Z)}-\theta ,\\ \eta (Z)& :=(g(0,Z),g(1,Z),m(Z)),{\eta }_0(Z):=(g_0(0,Z),g_0(1,Z),m_0(Z)),\end{array}$$

$$\begin{array}{rl}& \psi (W;\theta ,\eta )=\frac{D(Y-g(0,Z))}{m}-\frac{m(Z)(1-D)(Y-g(0,Z))}{(1-m(Z))m}-\theta \frac{D}{m},\\ \eta (Z)& :=(g(0,Z),g(1,Z),m(Z),m),{\eta }_0(Z)=(g_0(0,Z),g_0(1,Z),m_0(Z),\mathrm{E}[D]),\end{array}$$

$${\check{\theta }}_0(I_k,I_k^c),k=1,...,K,$$

$${\hat{\eta }}_0(I_k^c)={\left({\hat{g}}_0(0,Z;I_k^c),{\hat{g}}_0(1,Z;I_k^c),{\hat{m}}_0(Z;I_k^c),\frac{1}{N-n}\sum _{i\in I_k^c}{D}_i\right)}^{\prime },$$

$${\eta }_0(Z)=(g_0(0,Z),g_0(1,Z),m_0(Z),\text{E}[D]{)}^{\prime },$$

$$\frac{1}{n}\sum _{i\in I_k}\psi (W;\theta ,{\hat{\eta }}_0(I_k^c))=0,$$

$${\tilde{\theta }}_0=\frac{1}{K}\sum _{k=1}^K{\check{\theta }}_0(I_k,I_k^c).$$

$${\hat{\sigma }}^2=\frac{1}{N}\sum _{i=1}^N{\hat{\psi }}_i^2,$$

$${\text{CI}}_n:=[{\tilde{\theta }}_0\pm {\Phi }^{-1}(1-\alpha /2)\hat{\sigma }/\sqrt{N}].$$

$$\sqrt{s^g/n}\sqrt{s^m/n}\ll n^{-1/2}\iff s^g{s}^m\ll \sqrt{n}$$

$${\sigma }^{-1}\sqrt{N}({\tilde{\theta }}_0-{\theta }_0)⇝N(0,1),$$

$${\sigma }^2={\mathrm{E}}_P[{\psi }^2(W;{\theta }_0,{\eta }_0(Z))],$$

$$\underset{P\in \mathcal{P}}{\sup }|P\left({\theta }_0\in {\text{CI}}_n\right)-(1-\alpha )|\to 0.$$

$${\hat{\sigma }}^{\text{Mean}}=\sqrt{\frac{1}{S}\sum _{s=1}^S\left({\hat{\sigma }}_s^2+({\tilde{\theta }}_0^s-\frac{1}{S}\sum _{j=1}^S{\tilde{\theta }}_0^j{)}^2\right)},$$

$$\displaystyle\widehat{\sigma}^{\text{Median}}=\underset{}{\text{median}}\Big{\{}$$

$$\Vert f(W){\Vert }_{P,q}:={\left(\int |f(w){|}^{q}dP(w)\right)}^{1/q}.$$

$${\mathbb{G}}_{n,I}(\psi (W)):=\frac{1}{\sqrt{n}}\sum _{i\in I}\left(f(W_i)-\int f(w)dP(w)\right).$$

$${\mathbb{E}}_{n,I}(\psi (W)):=\frac{1}{n}\sum _{i\in I}f(W_i);{\mathbb{P}}_{n,I}(A):=\frac{1}{n}\sum _{i\in I}1(W_i\in A).$$

$$\begin{array}[]{ll}\psi(W;\theta,\eta):=g(1,Z)-g(0,Z)+\frac{D(Y-g(1,Z))}{m(Z)}-\frac{(1-D)(Y-g(0,Z))}{1-m(Z)}-\theta,\end{array}$$

$$\sqrt{n}({\check{\theta }}_k-{\theta }_0)={\mathbb{G}}_{n,I_k}\psi (W;{\theta }_0,{\hat{\eta }}_0(I_k^c))+\sqrt{n}\int \psi (w;{\theta }_0,{\hat{\eta }}_0(I_k^c))d{P}_n(w).$$

$$\int (\psi (w;{\theta }_0,{\hat{\eta }}_0(I_k^c))-\psi (w;{\theta }_0,{\eta }_0){)}^{2}d{P}_n(w)=o_{P_n}(1),$$

$$\sqrt{n}\int \psi (w;{\theta }_0,{\hat{\eta }}_0(I_k^c))d{P}_n(w)=o_{P_n}(1),$$

$${\hat{\sigma }}^2-{\sigma }^2=o_{P_n}(1),$$

$${\sigma }^{-2}=O(1),\text{and}{\sigma }^2=O(1).$$

$${\mathbb{G}}_{n,I_k}(\psi (W;{\theta }_0,{\hat{\eta }}_0(I_k^c))-\psi (W;{\theta }_0,{\eta }_0))=o_{P_n}(1),$$

$${\left({\sigma }^{-1}\sqrt{n}({\check{\theta }}_{0,k}-{\theta }_0)\right)}_{k=1}^K={\left({\sigma }^{-1}{\mathbb{G}}_{n,I_k}\psi (W;{\theta }_0,{\eta }_0)\right)}_{k=1}^K+o_{P_n}(1)⇝({\mathcal{N}}_k{)}_{k=1}^K,$$

$${\sigma }^{-1}\sqrt{nK}({\tilde{\theta }}_0-{\theta }_0)$$

Full text

Double/Debiased/Neyman Machine Learning of Treatment Effects

by

This is a conference version of \citeasnounCCDHM16. To appear in: American Economic Review 2017 (May), Papers & Proceedings.

Abstract. \citeasnounCCDHM16 provide a generic double/de-biased machine learning (DML) approach for obtaining valid inferential statements about focal parameters, using Neyman-orthogonal scores and cross-fitting, in settings where nuisance parameters are estimated using a new generation of nonparametric fitting methods for high-dimensional data, called machine learning methods. In this note, we illustrate the application of this method in the context of estimating average treatment effects (ATE) and average treatment effects on the treated (ATTE) using observational data. A more general discussion and references to the existing literature are available in \citeasnounCCDHM16.

**Key words: Neyman machine learning, orthogonalization, cross-fitting, double or de-biased machine learning, orthogonal score, efficient score, post-machine-learning and post-regularization inference, random forest, lasso, deep learning, neural nets, boosted trees, efficiency, optimality.

**

1. Scores for Average Treatment Effects

We consider estimation of ATE and ATTE under the unconfoundedness assumption of \citeasnounRR:prop. We consider the case where treatment effects are fully heterogeneous and the treatment variable, $D$, is binary, $D\in\{0,1\}$. We let $Y$ denote the outcome variable of interest and $Z$ denote a set of control variables. We then model random vector $(Y,D,Z)$ as

[TABLE]

Since $D$ is not additively separable, this model allows for very general heterogeneity in treatment effects. Common target parameters $\theta_{0}$ in this model are the ATE,

[TABLE]

and the ATTE,

[TABLE]

The confounding factors $Z$ affect the treatment variable $D$ via the propensity score, $m_{0}(Z):=\text{E}[D|Z],$ and the outcome variable via the function $g_{0}(D,Z)$. Both of these functions are unknown and potentially complicated, and we consider estimating these functions via the use of ML methods.

We proceed to set up moment conditions with scores that obey a type of orthogonality with respect to nuisance functions. Specifically, we make use of scores $\psi(W;\theta,\eta)$ that satisfy the identification condition

[TABLE]

and the Neyman orthogonality condition

[TABLE]

where $W=(Y,D,Z)$, $\theta_{0}$ is the parameter of interest, $\eta$ denotes nuisance functions with population value $\eta_{0}$, $\partial_{\eta}f\ |_{\eta=\eta_{0}}$ denote the derivative of $f$ with respect to $\eta$ (the Gateaux derivative operator).

Using moment conditions that satisfy (1.4) to construct estimators and inference procedures that are robust to small mistakes in nuisance parameters has a long history in statistics, following foundational work of \citeasnounNeyman59. Using moment conditions that satisfy (1.4) is also crucial to developing valid inference procedures for $\theta_{0}$ after using ML methods to produce estimators $\widehat{\eta}_{0}$ as discussed, e.g., in \citeasnounCHS:AnnRev. In practice, estimation of $\theta_{0}$ will be based on the empirical analog of (1.3) with $\eta_{0}$ replaced by $\widehat{\eta}_{0}$, and the Neyman orthogonality condition (1.4) ensures sufficient insensitivity to this replacement that high-quality inference for $\theta_{0}$ may be obtained. The second critical ingredient, that enables the use of wide array of modern ML estimators is data splitting, as discussed in the next section.

Neyman-orthogonal scores are readily available for both the ATE and ATTE – one can employ the doubly robust/efficient scores of \citeasnounrobins:dr and \citeasnounhahn:prop, which are automatically Neyman orthogonal. For estimating the ATE, we employ

[TABLE]

where $\eta(Z)$ is the nuisance parameter with true value denoted by $\eta_{0}(Z)$ consisting of $P$-square integrable functions, for $P$ defined in Assumption 2.1, mapping the support of $Z$ to $\mathbb{R}\times\mathbb{R}\times(\varepsilon,1-\varepsilon)$ where $\varepsilon>0$ is a constant. For estimation of ATTE, we use the score

[TABLE]

where again $\eta(Z)$ is the nuisance parameter with true value denoted by $\eta_{0}(Z)$ consisting of three $P$-square integrable functions, for $P$ defined in Assumption 2.1, mapping the support of $Z$ to $\mathbb{R}\times\mathbb{R}\times(\varepsilon,1-\varepsilon)$ and a constant $m\in(\varepsilon,1-\varepsilon)$. The respective scores for ATE and ATTE obey the identification condition (1.3) and the Neyman orthogonality property (1.4). Note that all semi-parametrically efficient scores share the orthogonality property (1.4), but not all orthogonal scores are efficient. In some problems, we may use inefficient orthogonal scores to have more robustness. Moreover, the use of efficient scores could be considerably refined using the targeted maximum likelihood approach of \citeasnounsrr:rejoinder and \citeasnounvdl:rubin in many contexts.

2. Algorithm and Result

We describe the DML estimator of $\theta_{0}$ using random sample $(W_{i})_{i=1}^{N}$. The algorithm makes use of a form of sample splitting, which we call cross-fitting. It builds on the ideas e.g. in \citeasnounssiv. The use of sample-splitting is a crucial ingredient to the approach that helps avoid overfitting which can easily result from the application of complex, flexible methods such as boosted linear and tree models, random forests, and various ensemble and hybrid ML methods. While the use of the Neyman orthogonality serves to reduce regularization and modeling biases of ML estimators $\widehat{\eta}_{0}$, the data splitting serves to eliminate the overfitting bias, which may easily occur, with the data scientist being even unaware of it. Thus the DML method employs the double debiasing.

Algorithm: Estimation using Orthogonal Scores by $K$-fold Cross-Fitting

Step 1. Let $K$ be a fixed integer. Form a $K$-fold random partition of $\{1,...,N\}$ by dividing it into equal parts $(I_{k})_{k=1}^{K}$ each of size $n:=N/K$, assuming that $N$ is a multiple of $K$. For each set $I_{k}$, let $I_{k}^{c}$ denote all observation indices that are not in $I_{k}$.

Step 2. Construct $K$ estimators

[TABLE]

that employ the machine learning estimators

[TABLE]

of the nuisance parameters

[TABLE]

and where each estimator $\check{\theta}_{0}(I_{k},I_{k}^{c})$ is defined as the root $\theta$ of

[TABLE]

for the score $\psi$ defined in (1.5) for the ATE and in (LABEL:ATTE-setup) for the ATTE.

Step 3. Average the $K$ estimators to obtain the final estimator:

[TABLE]

An approximate standard error for this estimator is $\widehat{\sigma}/\sqrt{N}$, where

[TABLE]

$\widehat{\psi}_{i}:=\psi(W_{i};\tilde{\theta}_{0},\widehat{\eta}_{0}(I^{c}_{k(i)}))$, and $k(i):=\{k\in\{1,...,K\}:i\in I_{k}\}$. An approximate $(1-\alpha)\times 100\%$ confidence interval is

[TABLE]

We now state a formal result that provides the asymptotic properties of $\tilde{\theta}_{0}$. Let $(\delta_{n})_{n=1}^{\infty}$ and $(\Delta_{n})_{n=1}^{\infty}$ be sequences of positive constants approaching 0. Let $c,\varepsilon,C$ and $q>4$ be fixed positive constants, and let $K$ be a fixed integer.

Assumption 2.1**.**

Let $\mathcal{P}$ be the set of probability distributions $P$ for $(Y,D,Z)$ such that (i) equations (1.1)-(1.2) hold, with $D\in\{0,1\}$, (ii) the following conditions on moments hold for all $d\in\{0,1\}$: $\|g(d,Z)\|_{P,q}\leqslant C$, $\|Y\|_{P,q}\leqslant C$, $P(\varepsilon\leqslant m_{0}(Z)\leqslant 1-\varepsilon)=1$, $P({\mathrm{E}}_{P}[\zeta^{2}\mid Z]\leqslant C)=1$, $\|\zeta\|_{P,2}\geqslant c$, and $\|\nu\|_{P,2}\geqslant c$ and (iii) the ML estimators of the nuisance parameters based upon a random subset $I^{c}_{k}$ of $\{1,...,N\}$ of size $N-n$, obey the following conditions for all $N\geqslant 2K$ and $d\in\{0,1\}$: $\|\widehat{g}_{0}(d,Z;I^{c}_{k})-g_{0}(d,Z)\|_{P,2}\cdot\|\widehat{m}_{0}(Z;I^{c}_{k})-m_{0}(Z)\|_{P,2}\leqslant\delta_{n}n^{-1/2}$, $\|\widehat{g}_{0}(d,Z;I^{c}_{k})-g_{0}(d,Z)\|_{P,2}+\|\widehat{m}_{0}(Z;I^{c}_{k})-m_{0}(Z)\|_{P,2}\leqslant\delta_{n}$, and $P(\varepsilon\leqslant\widehat{m}_{0}(Z;I^{c}_{k})\leqslant 1-\varepsilon)=1$, with ${\mathrm{P}}_{P}$-probability no less than $1-\Delta_{n}$.

The assumption on the rate of estimating the nuisance parameters is a non-primitive condition. These rates of convergence are available for most often used ML methods and are case-specific, so we do not restate conditions that are needed to reach these rates.

Comment 2.1** (Tighness of conditions).**

The conditions are fairly sharp (though somewhat simplified for presentation sake). The sharpness can be understood by examining the case where regression function $g_{0}$ and propensity function $m_{0}$ are sparse with sparsity indices $s^{g}\ll n$ and $s^{m}\ll n$, and estimators $\widehat{g}_{0}$ and $\widehat{m}_{0}$ having sparsity indices of orders $s^{g}$ and $s^{m}$ and covering to $g_{0}$ and $m_{0}$ at the rates $\sqrt{{s^{g}}/{n}}$ and $\sqrt{{s^{m}}/{n}}$, ignoring logs, based on $\ell_{1}$-penalized estimators. Then the rate conditions in the assumption require

[TABLE]

(ignoring logs) which is much weaker than the condition $(s^{g})^{2}+(s^{m})^{2}\ll n$ (ignoring logs) terms, required without sample splitting. For example, if the propensity score $m_{0}$ is very sparse, then the regression function is allowed to be quite dense, and vice versa. If the propensity score is known ($s^{m}=0$), then only consistency for $\widehat{g}_{0}$ is needed. Such comparisons also extend to the approximately sparse models. $\blacksquare$

Theorem 2.1**.**

Suppose that the ATE, $\theta_{0}=\mathrm{E}_{P}[g_{0}(1,Z)-g_{0}(0,Z)]$, is the target parameter and we use the estimator $\tilde{\theta}_{0}$ and other notations defined above. Alternatively, suppose that the ATTE, $\theta_{0}=\mathrm{E}_{P}[g_{0}(1,Z)-g_{0}(0,Z)\mid D=1]$, is the target parameter and we use the estimator $\tilde{\theta}_{0}$ and other notations above. Consider the set $\mathcal{P}$ of probability distributions $P$ defined in Assumption 2.1. Then, uniformly in $P\in\mathcal{P}$, the estimator $\tilde{\theta}_{0}$ concentrates around $\theta_{0}$ with the rate $1/\sqrt{N}$ and is approximately unbiased and normally distributed:

[TABLE]

and the result continues to hold if $\sigma^{2}$ is replaced by $\widehat{\sigma}^{2}$. Moreover, confidence regions based upon $\tilde{\theta}_{0}$ have uniform asymptotic validity:

[TABLE]

The scores $\psi$ are the efficient scores, so both estimators are asymptotically efficient, in the sense of reaching the semi-parametric efficiency bound of \citeasnounhahn:prop.

The proof follows from the application of Chebyshev’s inequality and the central limit theorem (and is given in the appendix for completeness).

3. Accounting for Uncertainty Due to Sample-Splitting

The method outlined in this note relies on subsampling to form auxiliary samples for estimating nuisance functions and main samples for estimating the parameter of interest. The specific sample partition has no impact on estimation results asymptotically but may be important in finite samples. Specifically, the dependence of the estimator on the particular split creates an additional source of variation. Incorporating a measure of this additional source of variation into estimated standard errors of parameters of interest may be important for quantifying the true uncertainty of the parameter estimates.

Hence we suggest making a slight modification to the asymptotically valid estimation procedure detailed in Section 2. Specifically, we propose repeating the main estimation procedure $S$ times, for a large number $S$, repartitioning the data in each replication $s=1,...,S$. Within each partition, we then obtain an estimate of the parameter of interest, $\tilde{\theta}_{0}^{s}$. Rather than report point estimates and interval estimates based on a single replication, we may then report estimates that incorporate information from the distribution of the individual estimates obtained from the $S$ different data partitions.

For point estimation, two natural quantities that could be reported are the sample average and the sample median of the estimates obtained across the $S$ replications, $\tilde{\theta}_{0}^{\text{Mean}}$ and $\tilde{\theta}_{0}^{\text{Median}}$. Both of these reduce the sensitivity of the estimate for $\theta_{0}$ to particular splits. $\tilde{\theta}_{0}^{\text{Mean}}$ could be strongly affected by any extreme point estimates obtained in the different random partitions of the data, and $\tilde{\theta}_{0}^{\text{Median}}$ is obviously much more robust. We note that asymptotically the specific random partition is irrelevant, and $\tilde{\theta}_{0}^{\text{Mean}}$ and $\tilde{\theta}_{0}^{\text{Median}}$ should be close to each other.

To quantify and incorporate the variation introduced by sample splitting, one might also compute standard errors that add an element to capture the spread of the estimates obtained across the $S$ different sets of partitions. For $\tilde{\theta}_{0}^{\text{Mean}}$, we propose adding an element that captures the spread of the estimated $\tilde{\theta}_{0}^{s}$ around $\tilde{\theta}_{0}^{\text{Mean}}$. Specifically, we suggest

[TABLE]

where $\widehat{\sigma}_{s}$ is defined as in Section 2. The second term in this formula takes into account the variation due to sample splitting which is added to a usual estimate of sampling uncertainty. Using this estimated standard error obviously results in more conservative inference than relying on the $\widehat{\sigma}_{s}$ alone. We adopt a similar formulation for $\tilde{\theta}_{0}^{\text{Median}}$. Specifically, we propose a median deviation defined as

[TABLE]

This standard error is more robust to outliers than $\widehat{\sigma}^{\text{Mean}}$.

Appendix A Practical Implementation and Empirical Examples

To illustrate the methods developed in this paper, we consider two empirical examples. In the first, we use the method outlined in the paper to estimate the effect of 401(k) eligibility on accumulated assets. In this example, the treatment variable is not randomly assigned and we aim to eliminate the potential biases due to the lack of random assignment by flexibly controlling for a rich set of variables. The second example reexamines the Pennsylvania Reemployment Bonus experiment which used a randomized control trial to investigate the incentive effect of unemployment insurance. Our goal in this supplement is to illustrate the use of our method and examine its empirical properties in two different settings: 1) an observational study where it is important to flexibly control for a large number of variables in order to overcome endogeneity, and 2) a randomized control trial where controlling for confounding factors is not needed for bias reduction but may produce more precise estimates.

Practical Details: Incorporating Uncertainty Induced by Sample Splitting

The results we report are based on repeated application of the method developed in the main paper as discussed in Section III. Specifically, we repeat the main estimation procedure 100 times repartitioning the data in each replication. We then report the average of the ATE estimates from the 100 random splits as the “Mean ATE,” and we report the median of the ATE estimates from the 100 splits as the “Median ATE.” We then report the measures of uncertainty that account for sampling variability and variability across the splits, $\widehat{\sigma}^{\text{Mean ATE}}$ and $\widehat{\sigma}^{\text{Median ATE}}$, for the “Mean ATE” and the “Median ATE” respectively.

The effect of 401(k) Eligibility on Net Financial Assets

The key problem in determining the effect of 401(k) eligibility is that working for a firm that offers access to a 401(k) plan is not randomly assigned. To overcome the lack of random assignment, we follow the strategy developed in \citeasnounpvw:94 and \citeasnounpvw:95. In these papers, the authors use data from the 1991 Survey of Income and Program Participation and argue that eligibility for enrolling in a 401(k) plan in this data can be taken as exogenous after conditioning on a few observables of which the most important for their argument is income. The basic idea of their argument is that, at least around the time 401(k)’s initially became available, people were unlikely to be basing their employment decisions on whether an employer offered a 401(k) but would instead focus on income and other aspects of the job. Following this argument, whether one is eligible for a 401(k) may then be taken as exogenous after appropriately conditioning on income and other control variables related to job choice.

A key component of the argument underlying the exogeneity of 401(k) eligibility is that eligibility may only be taken as exogenous after conditioning on income and other variables related to job choice that may correlate with whether a firm offers a 401(k). \citeasnounpvw:94 and \citeasnounpvw:95 and many subsequent papers adopt this argument but control only linearly for a small number of terms. One might wonder whether such specifications are able to adequately control for income and other related confounds. At the same time, the power to learn about treatment effects decreases as one allows more flexible models. The principled use of flexible machine learning tools offers one resolution to this tension. The results presented below thus complement previous results which rely on the assumption that confounding effects can adequately be controlled for by a small number of variables chosen ex ante by the researcher.

In the example in this paper, we use the same data as in \citeasnounCH401k. We use net financial assets - defined as the sum of IRA balances, 401(k) balances, checking accounts, U.S. saving bonds, other interest-earning accounts in banks and other financial institutions, other interest-earning assets (such as bonds held personally), stocks, and mutual funds less non-mortgage debt - as the outcome variable, $Y$, in our analysis. Our treatment variable, $D$, is an indicator for being eligible to enroll in a 401(k) plan. The vector of raw covariates, $Z$, consists of age, income, family size, years of education, a married indicator, a two-earner status indicator, a defined benefit pension status indicator, an IRA participation indicator, and a home ownership indicator.

In Table 1, we report estimates of the mean average treatment effect (Mean ATE) of 401(k) eligibility on net financial assets both in the partially linear model and allowing for heterogeneous treatment effects using the interactive model outlined in Section discussed in the main text. To reduce the disproportionate impact of extreme propensity score weights in the interactive model we trim the propensity scores which are close to the bounds, with the cutoff points of 0.01 and 0.99. We present two sets of results based on sample-splitting using a 2-fold cross-fitting and 5-fold cross-fitting.

We report results based on five simple methods for estimating the nuisance functions used in forming the orthogonal estimating equations. We consider three tree-based methods, labeled “Random Forest”, “Reg. Tree”, and “Boosting”, one $\ell_{1}$-penalization based method, labeled “Lasso”, and a neural network method, labeled “Neural Net”. For “Reg. Tree,” we fit a single CART tree to estimate each nuisance function with penalty parameter chosen by 10-fold cross-validation. The results in the “Random Forest” column are obtained by estimating each nuisance function with a random forest which averages over 1000 trees. The results in “Boosting” are obtained using boosted regression trees with regularization parameters chosen by 10-fold cross-validation. To estimate the nuisance functions using the neural networks, we use 8 hidden layers and a decay parameter of 0.01, and we set activation function as logistic for classification problems and as linear for regression problems.111We also experimented with “Deep Learning” methods from which we obtained similar results for some tuning parameters. However, we ran into stability and computational issues and chose not to report these results in the empirical section. “Lasso” estimates an $\ell_{1}$-penalized linear regression model using the data-driven penalty parameter selection rule developed in \citeasnounBellChenChernHans:nonGauss. For “Lasso”, we use a set of 275 potential control variables formed from the raw set of covariates and all second order terms, i.e. all squares and first-order interactions. For the remaining methods, we use the raw set of covariates as features.

We also consider two hybrid methods labeled “Ensemble” and “Best”. “Ensemble” optimally combines four of the machine learning methods listed above by estimating the nuisance functions as weighted averages of estimates from “Lasso,” “Boosting,” “Random Forest,” and “Neural Net”. The weights are restricted to sum to one and are chosen so that the weighted average of these methods gives the lowest average mean squared out-of-sample prediction error estimated using 5-fold cross-validation. The final column in Table 1 (“Best”) reports results that combines the methods in a different way. After obtaining estimates from the five simple methods and “Ensemble”, we select the best methods for estimating each nuisance functions based on the average out-of-sample prediction performance for the target variable associated to each nuisance function obtained from each of the previously described approaches. As a result, the reported estimate in the last column uses different machine learning methods to estimate different nuisance functions. Note that if a single method outperformed all the others in terms of prediction accuracy for all nuisance functions, the estimate in the “Best” column would be identical to the estimate reported under that method.

Turning to the results, it is first worth noting that the estimated ATE of 401(k) eligibility on net financial assets is $19,559 (not reported) with an estimated standard error of 1413 when no control variables are used. Of course, this number is not a valid estimate of the causal effect of 401(k) eligibility on financial assets if there are neglected confounding variables as suggested by \citeasnounpvw:94 and \citeasnounpvw:95. When we turn to the estimates that flexibly account for confounding reported in Table 1, we see that they are substantially attenuated relative to this baseline that does not account for confounding, suggesting much smaller causal effects of 401(k) eligiblity on financial asset holdings. It is interesting and reassuring that the results obtained from the different flexible methods are broadly consistent with each other. This similarity is consistent with the theory that suggests that results obtained through the use of orthogonal estimating equations and any sensible method of estimating the necessary nuisance functions should be similar. Finally, it is interesting that these results are also broadly consistent with those reported in the original work of \citeasnounpvw:94 and \citeasnounpvw:95 which used a simple intuitively motivated functional form, suggesting that this intuitive choice was sufficiently flexible to capture much of the confounding variation in this example.

There are other interesting observations that can provide useful insights into understanding the finite sample properties of the “double ML” estimation method. First, the standard errors of the estimates obtained using 5-fold cross-fitting are considerably lower than those obtained from 2-fold cross-fitting for all methods. This fact suggests that having more observations in the auxiliary sample may be desirable. Specifically, 5-fold cross-fitting estimates uses more observations to learn the nuisance functions than 2-fold cross-fitting and thus likely learns them more precisely. This increase in precision in learning the nuisance functions may then translate into more precisely estimated parameters of interest. While intuitive, we note that this statement does not seem to be generalizable in that there does not appear to be a general relationship between the number of folds in cross-fitting and the precision of the estimate of the parameter of interest. Second, we also see that the standard errors of the Lasso estimates are noticeably larger than the standard errors coming from the other machine learning methods. We believe that this is due to the fact that the out-of-sample prediction errors from a linear model tend to be larger when there is a need to extrapolate. In our framework, if the main sample includes observations that are outside of the range of the observations in the auxiliary sample, the model has to extrapolate to those observations. The fact that the standard errors are lower in 5-fold cross-fitting than in 2-fold cross-fitting for the “Lasso” estimations also supports this hypothesis, because the higher number of observations in the auxiliary sample reduces the degree of extrapolation.

In Table 2, we report the median ATE estimation results for the same models to check the robustness of our results to the outliers due to sample splitting. We see that both the coefficients and standard errors are similar to the “Mean ATE” estimates and standard errors. The similarity between the “Mean ATE” and “Median ATE” suggests that the distribution of the ATE across different splits is approximately symmetric and relatively thin-tailed.

The effect of Unemployment Insurance Bonus on Unemployment Duration

As a further example, we re-analyze the Pennsylvania Reemployment Bonus experiment which was conducted by the US Department of Labor in the 1980s to test the incentive effects of alternative compensation schemes for unemployment insurance (UI). This experiment has been previously studied by \citeasnounB:PennSeq and \citeasnounBK:PennQuan. In these experiments, UI claimants were randomly assigned either to a control group or one of five treatment groups.222There are six treatment groups in the experiments. Following \citeasnounB:PennSeq. we merge the groups 4 and 6. In the control group the standard rules of the UI applied. Individuals in the treatment groups were offered a cash bonus if they found a job within some pre-specified period of time (qualification period), provided that the job was retained for a specified duration. The treatments differed in the level of the bonus, the length of the qualification period, and whether the bonus was declining over time in the qualification period; see \citeasnounBK:PennQuan for further details on data.

In our empirical example, we focus only on the most generous compensation scheme, treatment 4, and drop all individuals who received other treatments. In this treatment, the bonus amount is high and qualification period is long compared to other treatments, and claimants are eligible to enroll in a workshop. Our treatment variable, D, is an indicator variable for the treatment 4, and the outcome variable, Y, is the log of duration of unemployment for the UI claimants. The vector of covariates, Z, consists of age group dummies, gender, race, number of dependents, quarter of the experiment, location within the state, existence of recall expectations, and type of occupation.

We report estimates of the ATE on unemployment duration both in the partially linear model and in the interactive model. We again consider the same methods with the same tuning choices, with one exception, for estimating the nuisance functions as in the previous example and so do not repeat details for brevity. The one exception is that we implement neural networks with 2 hidden layers and a decay parameter of 0.02 in this example which yields better prediction performance. In “Lasso” estimation, we use a set of 96 control variables formed by taking nonlinear functions and interactions of the raw set of covariates. For the remaining approaches, we use only the 14 raw control variables listed above.

Table 3 presents estimates of the “Mean ATE” on unemployment duration using the partially linear model and interactive model in panel A and B, respectively. To reduce the disproportionate impact of extreme propensity score weights in the interactive model, we trim the propensity scores at 0.01 and 0.99 as in the previous example. For both the partially linear model and the interactive model, we report estimates obtained using 2-fold cross-fitting and 5-fold cross-fitting.

The estimation results are consistent with the findings of previous studies which have analyzed the Pennsylvania Bonus Experiment. The ATE on unemployment duration is negative and significant across all estimation methods at the 5% level with the exception of the estimate of the ATE obtained from the interactive model using random forests, which is significant at the 10% level. When looking at standard errors it is useful to remember that they include both sampling variation and variation due to random sample splitting. It is reassuring to see that the variation due to sample splitting does not change the conclusion. It is also interesting to see that, similar to the result in the first empirical example, the “Mean ATE” estimates are broadly similar across different estimation models. Finally in Table 4 we report the “Median ATE” estimates. The median estimates are close to the mean estimates, giving further evidence for the stability of estimation across different random splits.

In conclusion, we want to emphasize some important observations that can be drawn from these empirical examples. First, for both examples the choice of the machine learning method in estimating nuisance functions does not substantively change the conclusion, and we obtained broadly consistent results regardless of which method we employ. Second, the similarity between the median and mean estimates suggests that the results are robust to the particular sample split used in estimation in these examples.

Appendix B Proofs

Notation

The symbols ${\mathrm{P}}$ and ${\mathrm{E}}$ denote probability and expectation operators with respect to a generic probability measure. If we need to signify the dependence on a probability measure $P$, we use $P$ as a subscript in ${\mathrm{P}}_{P}$ and ${\mathrm{E}}_{P}$. Note also that we use capital letters such as $W$ to denote random elements and use the corresponding lower case letters such as $w$ to denote fixed values that these random elements can take in the set $\mathcal{W}$. In what follows, we use $\|\cdot\|_{P,q}$ to denote the $L^{q}(P)$ norm; for example, for measurable $f:\mathcal{W}\to\mathbb{R}$, we denote

[TABLE]

Define the empirical process $\mathbb{G}_{n,I}(\psi(W))$ as a linear operator acting on measurable functions $\psi:\mathcal{W}\to\mathbb{R}$ such that $\|\psi(W)\|_{P,2}<\infty$ via,

[TABLE]

Analogously, we define the empirical expectation and probability as:

[TABLE]

Proof of Theorem 1

We will demonstrate the result for the case of ATE estimator, which uses the score:

[TABLE]

and the result for ATTE follows similarly. Choose any sequence $\{P_{n}\}\in\mathcal{P}$.

Step 1: (Main Step). Letting $\check{\theta}_{0,k}=\check{\theta}_{0}(I_{k},I_{k}^{c})$, write

[TABLE]

Steps 2 – 5 below demonstrate that for each $k=1,...,K$,

[TABLE]

The equations (B.1) and (B.2) are the minimal conditions needed on the estimators of the nuisance parameters, and could be used to replace the more primitive conditions stated in the text.

Assertion (B.1) implies that

[TABLE]

since the quantity in the display converges in probability conditionally on the data $(W_{i})_{i\in I^{c}_{k}}$ by (B.1) and Chebyshev inequality, which in turn implies the unconditional convergence in probability, as noted in the following simple lemma.

Lemma B.1**.**

Let $\{X_{m}\}$ and $\{Y_{m}\}$ be sequences of random vectors. If for any $\epsilon>0$, ${\mathrm{P}}(\|X_{m}\|>\epsilon\mid Y_{m})\to_{{\mathrm{P}}}0$, then ${\mathrm{P}}(\|X_{m}\|>\epsilon)\to 0$. In particular, this occurs if ${\mathrm{E}}[\|X_{m}\|^{q}\mid Y_{m}]\to_{{\mathrm{P}}}0$ for some $q\geqslant 1$, by Markov’s inequality.

Proof. For any $\epsilon>0$ ${\mathrm{P}}(\|X_{m}\|>\epsilon)\leqslant{\mathrm{E}}[{\mathrm{P}}(\|X_{m}\|>\epsilon\mid Y_{m})]\to 0$, since the sequence $\{{\mathrm{P}}(\|X_{m}\|>\epsilon\mid Y_{m})\}$ is uniformly integrable. $\blacksquare$

Using independence of data blocks $(W_{i})_{i\in I_{k}}$, $k=1,...,K$, the application of the Lindeberg-Feller theorem and the Cramer-Wold device, we conclude that

[TABLE]

where $(\mathcal{N}_{k})_{k=1}^{K}$ is a Gaussian vector with independent $N(0,1)$ coordinates. Therefore,

[TABLE]

where the last line uses the sum-stability of the normal distribution. Moreover, the result continues to hold if $\sigma$ is replaced by $\widehat{\sigma}$ in view of (B.3) and (B.4).

The above claim implies that $\mathrm{CI}_{n}=[\tilde{\theta}_{0}\pm\Phi^{-1}(1-\alpha/2)\widehat{\sigma}/\sqrt{N}]$ obeys

[TABLE]

In addition, the last two claims hold under any sequence $\{P_{n}\}\in\mathcal{P}$, which implies that these claims hold uniformly in $P\in\mathcal{P}$. Indeed, for example, choose $\{P_{n}\}$ such that, for some $\epsilon_{n}\to 0$

[TABLE]

The right-hand side converges to zero, which implies the uniform convergence. It remains to prove claims (B.1) – (B.4).

Step 2: This step demonstrates assertion (B.1). Observe that for some constant $C_{\varepsilon}$ that depends only on $\varepsilon$ and $\mathcal{P}$,

[TABLE]

where

[TABLE]

We bound $\mathcal{I}_{1}$, $\mathcal{I}_{2}$, and $\mathcal{I}_{3}$ in turn. First, $P_{n}(\mathcal{I}_{1}>\delta_{n})\leqslant\Delta_{n}\to 0$ by Assumption 2.1, and so $\mathcal{I}_{1}=o_{P_{n}}(1)$. Also, on the event that

[TABLE]

and

[TABLE]

which happens with ${\mathrm{P}}_{P_{n}}$-probability at least $1-\Delta_{n}$ by Assumption 2.1,

[TABLE]

where the first inequality follows from (B.5) and Assumption 2.1, the second from the facts that $D\in\{0,1\}$ and for $D=1$, $Y=g_{0}(1,Z)+\zeta$, the third from the triangle inequality, the fourth from the facts that $P_{n}(m_{0}(Z)\leqslant 1)=1$ and $P_{n}({\mathrm{E}}_{P_{n}}[\zeta^{2}\mid Z]\leqslant C)=1$, which are imposed in Assumption 2.1, the fifth from (B.6), and the last assertion follows since $\delta_{n}\to 0$. Hence, $\mathcal{I}_{2}=o_{P_{n}}(1)$. In addition, the same argument shows that $\mathcal{I}_{3}=o_{P_{n}}(1)$, and so the claim of this step follows.

Step 3: This step demonstrates the assertion (B.2). Observe that since $\theta_{0}={\mathrm{E}}_{P_{n}}[g_{0}(1,Z)-g_{0}(0,Z)]$, the left-hand side of (B.2) is equal to

[TABLE]

But on the event that

[TABLE]

and

[TABLE]

which happens with ${\mathrm{P}}_{P_{n}}$-probability at least $1-\Delta_{n}$ by Assumption 2.1, we have using the Cauchy-Schwarz inequality that

[TABLE]

which gives the claim of this step.

Step 4: This step demonstrates the assertion (B.3). We can write

[TABLE]

We claim that for each $k=1,...,K$,

[TABLE]

The latter property holds by the Chebyshev Inequality. Further, letting $I$ denote a generic $I_{k}$, the relation $a^{2}-b^{2}=(a-b)(a+b)$, and the Cauchy-Schwarz and triangle inequalities yield

[TABLE]

where

[TABLE]

Hence, given that $\|\psi(W;\theta_{0},\eta_{0})\|^{2}_{\mathbb{P}_{n,I},2}-\sigma^{2}=o_{P_{n}}(1)$ as noted above, and $\sigma^{2}$ is bounded above by (B.4), the claim follows as long as we can show that

[TABLE]

To show (B.7), we have

[TABLE]

and $\|\tilde{\theta}_{0}-\theta_{0}\|=o_{P_{n}}(1)$ by the arguments in Step 1. Also,

[TABLE]

by Jensen’s inequality and Step 2. Hence, (B.7) holds, which gives the claim of this step.

Step 5: This step demonstrates the assertion (B.4). Note that the bound $\sigma^{2}=O(1)$ holds trivially from Assumption 2.1. Hence, it suffices to show that $\sigma^{-2}=O(1)$. To do so, observe that using Assumption 2.1 we have

[TABLE]

where the last line follows since $P_{n}(m(Z)\leqslant 1)=P_{n}(1-m(Z)\leqslant 1)=1$. This gives the claim of this step and completes the proof of the theorem. $\blacksquare$