Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
77,223 characters · 20 sections · 48 citation commands
Doubly-Robust Inference for Conditional Average Treatment Effects with High-Dimensional Controls
Consider a potential outcomes framework RubinPotentialOutcomes, Rubin_1978_AnnalsStat where an observed outcome \(Y \in \SR\) and treatment \(D \in \{0,1\}\) are related to two latent potential outcomes \(Y_1, Y_0\in\SR\) via \(Y = DY_1 + (1-D)Y_0\). To account for unobserved confounding factors a common strategy is to assume the researcher has access to a vector of covariates, \(Z = (Z_1, X) \in \calZ_1 \times \calX \subseteq \SR^{d_z-d_x,d_x}\), such that the potential outcomes are independent of the treatment decision after conditioning on the observed covariates, \((Y_1,Y_0)\perp D | Z\). In this setting, we are interested in estimation of and inference on the conditional average treatment effect (CATE):
Estimation of the CATE generally requires first fitting propensity score and/or outcome regression models. When the number of control variables \(Z\) is large (\(d_z \gg n\)), these first stage models must be estimated using regularized methods which converge slower than the nonparametric rate and typically rely on the correctness of parametric specifications for consistency.\footnote{Recent works by BK-2019-neural-nets,Schmidt-Hieber-2020-neural-networks provide some limited nonparametric results in high-dimensional settings using deep neural networks.}
Fortunately, so long as both models are correctly specified, one can obtain a nonparametric-rate consistent estimator and valid inference procedure for the CATE by using the popular augmented inverse propensity weighted (aIPW) signal SC-2020,fan2022estimation. This is because the aIPW signal obeys an orthogonality condition at the true nuisance model values that limits the first stage estimation error passed on to the second stage estimator. Moreover, estimators based on the aIPW signal are doubly-robust; consistency of the resulting second-stage estimators requires correct specification of only one of the first stage propensity score or outcome regression models. However inference based on these estimators is not doubly-robust. Under misspecification the aIPW signal orthogonality fails and resulting testing procedures and confidence intervals are rendered invalid.
This paper proposes a doubly-robust estimator and inference procedure for the conditional average treatment effect when the number of control variables \(d_z\) is potentially much larger than the sample size \(n\). The dimensionality of the conditioning variable, \(d_x\), remains fixed in our analysis. Our approach is based on Tan-2018 wherein doubly-robust inference is developed for the average treatment effect. Following SC-2020 we take a series approach to estimating the CATE, using a quasi-projection of the aIPW signal onto a growing set of basis functions. By assuming a logistic form for the propensity score model and a linear form for the outcome regression model, we construct novel \(\ell_1\)-regularized first-stage estimating equations to recover a partial orthogonality of the aIPW signal at the limiting values of the first stage estimators. This restricted orthogonality is enough to achieve doubly robust pointwise and uniform inference; pointwise and uniform confidence intervals centered at the second-stage estimator are valid even if one of the logistic or linear functional forms is misspecified.
To achieve doubly-robust inference at all points in the support of the conditioning variable, we must obtain this restricted orthogonality for each basis term in the series approximation. This is accomplished by employing distinct first-stage estimating equations for each basis term used in the second-stage series approximation. This results in the number of first-stage estimators growing with the number of basis terms. These estimators converge uniformly to limiting values under standard conditions in high-dimensional analysis. Improving on prior work in doubly-robust inference, our \(\ell_1\) regularized first-stage estimation incorporates a data-dependent penalty parameter based on the work of CS-2021. This allows practical implementation of our proposed estimation procedure with minimal knowledge of the underlying data generating process.
The use of multiple pairs of nuisance parameter estimates limits our ability to straightforwardly apply existing nonparametric results for series estimators Newey-1997,BCCK-2015. Under modified conditions, we analyze the asymptotic properties of our second-stage series estimator to re-derive pointwise and uniform inference results. These modified conditions are in general slightly stronger than those of BCCK-2015, though in certain special cases collapse exactly to the conditions of BCCK-2015.
\paragraph{Prior Literature.} CCDDHNR-2018 analyze the general problem of estimating finite dimensional target parameters in the presence of potentially high dimensional nuisance functions. Using score functions that are Neyman-orthogonal with respect to nuisance parameters they show that it is possible to obtain target parameter estimates that are \(\sqrt{n}\)-consistent and asymptotically normal so long as the nuisance parameters are consistent at rate \(n^{-1/4}\), a condition satisfied by many machine learning-based estimators. SC-2020 take advantage of new results for series estimation in BCCK-2015 and consider series estimation of functional target parameters after high-dimensional nuisance estimation.\footnote{fan2022estimation provides a similar analysis using a second stage kernel estimator.}
In the same setting as this paper, Tan-2018 considers estimation of the average treatment effect. After assuming a logistic form for the propensity score and a linear form for the outcome regression, Tan-2018 proposes \(\ell_1\)-regularized first-stage estimators that allow for partial control of the derivative of the aIPW signal away from true nuisance values and thus allow for doubly-robust inference. SRR-2019 extends the analysis of Tan-2018 to consider doubly-robust inference for a larger class of finite dimensional target parameters with bilinear influence functions. tanCATE provide doubly-robust inference procedures for covariate-specific treatment effects with discrete conditioning variables; their results depend on exact representation assumptions that are unlikely to hold with continuous covariates. Moreover, no uniform inference procedures are described.
CS-2021 propose a data-driven “bootstrap after cross-validation” approach to penalty parameter selection that is modified for and implemented in our setting. This work is related to other work on the lasso tibshirani1996regression,BRT-2009,BC-2013,CLC-2021-CVLasso and \(\ell_1\)-regularized M-estimation in high dimensional settings vanDerGreer2016,Tan-2017.
\paragraph{Paper Structure.} This paper proceeds as follows. Section (ref) defines the problem and introduces our methods for estimation and inference. Section (ref) provides intuition for how the first stage estimation procedure allows for doubly-robust estimation and inference on the CATE as well as formally establishes the necessary first stage convergence. Section (ref) presents the main results: valid pointwise and uniform inference for the second-stage series estimator if either the first-stage logistic propensity score model or linear outcome regression model is correctly specified. Section (ref) ties up a technical detail. Section (ref) provides evidence from a simulation study while Section (ref) applies our proposed estimator to examine the effect of maternal smoking on infant birth weight. Section (ref) concludes. Proofs of main results are deferred to the Appendix.
\paragraph{Notation.} For any measure \(F\) and any function \(f\), define the \(L^2\) norm, \(\|f\|_{F, 2} = (\E_{F}[f^2])^{1/2}\) and the \(L^\infty\) norm \(\|f\|_{F, \infty} = \esssup_F |f| \). For any vector in \(\SR^p\) let \(\|\cdot\|_p\) for \(p \in [1,\infty]\) denote the \(\ell_p\) norm, \(\|a\|_p = (\sum_{l=1}^p a_l^p)^{1/p}\) and \(\|a\|_\infty = \max_{1\leq l\leq \infty}|a_l|\). If the subscript is unspecified, we are using the \(\ell_2\) norm. For two vectors \(a,b\in\SR^p\), let \(a\circ b = (a_ib_i)_{i=1}^p\) denote the Hadamard (element-wise) product. We adopt the convention that for \(a\in\SR^p\) and \(c\in\SR\), \(a + c = (a_i + c)_{i=1}^p\). For a matrix \(A\in\SR^{m\times n}\) let \(\|A\| = \max_{\|v\|_{\ell_2} \leq 1}\|Av\|_{\ell_2}\) denote the operator norm and \(\|A\|_\infty = \sup_{{1\leq r\leq m,1\leq s \leq n}} |A_{rs}|\). For any real valued function \(f\) let \(\E_n[f(X)] = \frac{1}{n} \sum_{i=1}^n f(X_i)\) denote the empirical expectation and \( \mathbb{G}_n[f(X)] = \frac{1}{\sqrt{n}}\sum_{i=1}^n (f(X_i) - \E[X_i])\) denote the empirical process. For two sequences of random variables \(\{a_n\}_{\SN}\) and \(\{b_n\}_{\SN}\), we say \(a_n \lesssim_P b_n\) or \(a_n = O_p(b_n)\) if \(a_n/b_n\) is bounded in probability and say \(a_n = o_p(b_n)\) if \(a_n/b_n \to_p 0\).
Below, we formally define the setting and identification strategy that we consider. We then introduce our doubly-robust estimator and inference procedure. The parameter of interest is the conditional average treatment effect: \(\E[Y_1 - Y_0\mid X=x]\). However, for this paper we largely focus on estimation and inference for the conditional average counterfactual outcome:
Doubly-robust estimation and inference on the other conditional counterfactual outcome, \(\E[Y_0\,|\,X=x]\), follows a similar procedure and is described in (ref). The procedures can be combined for doubly-robust estimation and inference for the CATE.
We assume that the researcher observes i.i.d data and that conditioning on \(Z\) is sufficient to control for all confounding factors affecting both the treatment decision \(D\) and the potential outcomes, \(Y_1\) and \(Y_0\). Our analysis allows the dimensionality of the controls, \(Z = (Z_1,X)\), to grow much faster than sample size \((d_z \gg n)\), while assuming the dimensionality of the conditioning variables, \(X\), remains fixed \((d_x \ll n)\).
To obtain doubly-robust estimation and inference we use the augmented inverse propensity weighted (aIPW) signal,
which is a function of a fitted propensity score model, \(\pi(Z),\) and a fitted outcome regression model, \(m(Z)\), whose true values are given \(\pi^\star(Z) := \E[D\mid Z]\) and \(m^\star(Z) := \E[Y\mid D= 1, Z]\). Under (ref), the aIPW signal \(Y(\cdot,\cdot)\) provides doubly-robust identification of \(g_0(x)\). That is, for integrable \(\pi\neq \pi^\star\) and \(m\neq m^\star\),
We use a series approach to estimate \(g_0(x)\), taking a quasi-projection of the aIPW signal onto a growing set of \(k\) weakly positive basis terms:
The basis terms are required to be weakly positive as they are used as weights within the convex first-stage estimators estimating equations.\footnotemark Examples of weakly positive basis functions are B-splines or shifted polynomial series terms. To ensure that the basis terms are well behaved, we make assumptions on \(\xi_{k,\infty} := \sup_{x\in\calX}\|p^k(x)\|_\infty\), \(\xi_{k,2} := \sup_{x\in\calX}\|p^k(x)\|_2\), and the eigenvalues of the design matrix \(Q := \E[p^k(x)p^k(x)']\).
For each basis term \(p_j(x), j = 1,\dots, k\), we estimate a separate propensity score model, \(\widehat\pi_j(Z)\), and outcome regression model, \(\widehat m_j(Z)\). Under standard moment and sparsity conditions, these converge uniformly over \(j=1,\dots,k\) to limiting values \(\bar\pi_j(Z)\) and \(\bar m_j(Z)\). If the propensity score model and outcome regression models are correctly specified these limiting values coincide with the true values \(\pi^\star(Z)\) and \(m^\star(Z)\). However, in general the limiting and true values may differ. The double robustness of the aIPW signal allows for identification of the CATE even if only one of the nuisance models is correctly specified. If either \(\bar\pi_j = \pi^\star\) or \(\bar m_j = m^\star\), we can write for all \(j=1,\dots,k\):
where \(g_0(x)\) is the conditional counterfactual outcome (ref), \(g_k(x) := p^k(x)'\beta^k\) is the projection of \(g_0(x)\) onto the first \(k\) basis terms, and \(r_k (x) := g_0(x) - g_k (x) \) denotes the approximation error from this projection. Note the separate error terms for each \(j=1,\dots,k\) in (ref), which are collected together in the vector \(\eps^k := (\eps_1,\dots,\eps_k)\). As long as one of the first-stage models is correctly specified, the least squares parameter \( \beta^k \) governing the projection in \(g_k(x)\) can be identified by the projection of the aIPW signal onto the basis terms \(p^k(x)\):
\footnotetext{(ref) provides a slightly modified method of constructing our doubly-robust estimator and inference procedure that does not require the first stage weights to directly be the second stage basis terms. This may be useful in case the researcher wants to use a second stage basis that cannot be transformed to be weakly positive.}
We assume a logistic regression form for the propensity score model and a linear form for the outcome regression model:
For each \(j=1,\dots,k,\) the parameters of (ref), \(\gamma,\alpha \in \SR^{d_z} ,\) are estimated by
The penalty parameters \(\lambda_{\gamma, j}\) and \(\alpha_{\gamma, j}\) are chosen via a data dependent technique described below. These first stage estimating equations are designed so that their first order conditions directly limit the bias passed on to the second-stage series estimator, as is described in (ref). Under standard assumptions the parameter estimators \(\widehat\gamma_j, \widehat\alpha_j\) will converge uniformly over \(j=1,\dots,k\) to population minimizers
which we assume are sufficiently sparse. Our first stage estimators are then \(\widehat\pi_j(Z) := \pi(Z;\widehat\gamma_j)\) and \(\widehat m_j(Z) := m(Z;\widehat\alpha_j)\) with limiting values \(\bar\pi_j(Z) := \pi(Z;\bar\gamma_j)\) and \(\bar m_j(Z) := m(Z;\bar\alpha_j)\), respectively.
Our second stage estimator is then \(\widehat g(x) := p^k(x)'\widehat\beta^k\) where \(\widehat\beta^k\) is an estimate of the population projection parameter, \(\beta^k\), obtained by combining all \(k\) pairs of first stage estimators according to
and \(\widehat Q := \E_n[p^k(X)p^k(X)']\). We estimate the variance of \(\widehat g(x)\) using \(\widehat\sigma(x) := \|\widehat\Omega^{1/2}p^k(x)\|/\sqrt{n}\) for
where \(\circ\) represents the Hadamard product and \(\widehat\eps^k := (\widehat\eps_1,\dots,\widehat\eps_k)\); \(\widehat\eps_j := Y(\widehat\pi_j, \widehat m_j) - \widehat g(x)\), \(j=1,...,k\).
Inference is based on the \(100(1-\eta)\%\) confidence bands
For pointwise inference, the critical value \(c^\star(1-\eta/2)\) is taken as the \((1-\eta/2)\) quantile of a standard normal distribution. For uniform inference \(c^\star(1-\eta/2)\) is taken \[ c^\star(1-\eta/2) := (1-\eta/2)\text{-quantile of }\sup_{x\in\calX} \left|\frac{p^k(x)\widehat\Omega^{1/2}}{\widehat\sigma(x)}N_k^b \right| \] where \(N_k^b\) is a bootstrap draw from \(N(0,I_k)\). (ref) show that, under standard sparsity and moment conditions, these pointwise and uniform inference procedures remain valid even under misspecification of either first-stage model.
To select the penalty parameters \(\lambda_{\gamma,j}\) and \(\lambda_{\alpha,j}\) in (ref)-(ref) we propose a data driven two-step procedure based on the work of CS-2021. For each \(j= 0,1\dots,k,\) we start with pilot penalty parameters given by
for some constants \(c_{\gamma, j}, c_{\alpha, j}\) selected from the interval \([\underline c_n, \bar c_n]\) with \(\underline c_n > 0\). In practice, the researcher has a fair bit of flexibility in choosing these constants. The optimal choice of these constants may depend on the underlying data generating process. We recommend using cross validation to pick these constants from a fixed-cardinality set of possible values. In line with (ref)(vi), the values in the set should be chosen to be on the order of the maximum value of \(\|p^k(X_i)\|_\infty\) observed in the data.
Using \(\lambda^{\text{\tiny pilot}}_{\gamma, j}\) and \(\lambda^{\text{\tiny pilot}}_{\alpha,j}\) in lieu of \(\lambda_{\gamma,j}\) and \(\lambda_{\alpha,j}\) in (ref)-(ref) we generate pilot estimators \(\widehat\gamma^{\text{\tiny pilot}}_j\) and \(\widehat\alpha^{\text{\tiny pilot}}_j\). These pilot estimators are used to generate plug in estimators \(\widehat U_{\gamma, j}\) and \(\widehat U_{\alpha, j}\) of the residuals
We then use a multiplier bootstrap procedure to select our final penalty parameters \(\lambda_{\gamma, j}\) and \(\lambda_{\alpha, j}\).
where \(e_1,\dots,e_n\) are independent standard normal random variables generated independently of the data \(\{Y_i,D_i,X_i\}_{i=1}^n\) and \(c_0 > 1\) is a fixed constant.\footnote{The constant \(c_0\) can be different for the propensity score and outcome regression models and can also vary for each \(j = 1,\dots,k\). All that matters is that each constant satisfies the requirements of (ref). This complicates notation, however.} In line with other work we find \(c_0 = 1.1\) works well in simulations. So long as our residual estimates converge in empirical mean square to limiting values, the choice of penalty parameter in (ref) will ensure that the penalty parameter dominates the noise with high probability. This allows for consistent variable selection and coefficient estimation.
For computational reasons, the researcher may not want to implement the bootstrap penalty parameter procedure. If this is the case, we note that the pilot penalty parameters of (ref) can be used directly after the constants \(c_{\gamma, j}\) and \(c_{\alpha,j}\) can be selected via cross validation from a growing set \(\Lambda_n \subseteq [\underline c_n, \bar c_n]\) under modified conditions. (ref) provides details for this implementation as well as formally shows the modified conditions needed.
We begin with a main technical lemma which provides a bound on rate at which first stage estimation error is passed on to the second stage CATE and variance estimators. This bound is comparable to others seen in the inference after model-selection literature BCH-2013,Tan-2018 and is achieved under standard conditions in the \(\ell_1\)-regularized estimation literature BRT-2009,Bulmann-VanDeGeer-2011,BC-2013,CS-2021. However, this bound is achieved at the limiting values of the propensity score and outcome regression models which may differ from the true values \(\pi^\star\) and \(m^\star\) under misspecification.
The potential misspecification of the first stage models which means we cannot directly apply orthogonality of the aIPW signal, discussed below, to show that the effect of first stage estimation error on the second stage is negligible. Instead, we use the first order conditions for \(\widehat\gamma_j\) and \(\widehat\alpha_j\) to directly control this quantity. After presenting the lemma (ref) provides some intuition for how this is done. Controlling the rate at which first stage estimation error is passed on to the second stage estimator even at points away from the true values \(\pi^\star\) and \(m^\star\) is key for obtaining doubly-robust inference for the CATE.
To show uniform convergence of the first stage estimators and thus uniform control of the bias passed on from the first stage estimation to the second stage estimator we rely on the following assumption:
The first part of (ref) assumes that the regressors are bounded while the second assumes that tail behavior of the outcome regression errors are uniformly thin. Both of these can be relaxed somewhat with sufficient moment conditions on the tail behavior of the controls and errors. We should note that compactness of \(\calX\) is generally required by nonparametric estimators. The third part of the assumption bounds all limiting propensity scores \(\bar\pi_j(Z)\) away from zero uniformly. The fourth assumption is an empirical compatibility condition on the weighted first-stage design matrix. It is slightly weaker than the restricted eigenvalue conditions often assumed in the literature BRT-2009,BCCK-2012. The penultimate condition is an identifiability constraint that limits the moments of the noise and bounds it away from zero uniformly over all estimation procedures. Many of the constants in (ref) are assumed to be fixed across all \(j\). This is mainly to simplify the exposition of the results below and in practice all constants can be allowed to grow slowly with \(k\). However, the growth rate of these terms affects the required first-stage sparsity.
The last condition is required for the validity of the bootstrap penalty parameter selection procedure and is comparable to the requirements needed for the bootstrap after cross validation technique described by CS-2021. The main difference is the additional assumption on the growth rate of the basis functions, \(\xi_{k,\infty}\) which is to ensure uniform stability of the estimation procedures (ref)-(ref) as well as some assumptions on the order of the constants \(c_{\gamma,j}\) and \(c_{\alpha,j}\) in (ref).
(ref) provides a tight bound on the first-stage estimation error passed on to the second stage estimator even when the first-stage estimators converge to values that are not the true propensity score or outcome regression. In particular notice that under the (familiar) sparsity bound \(s_k\xi_{k,\infty}^2k^{1/2}\ln^2(d_z)/\sqrt{n}\to 0\), any linear combination of the means in both (ref) and (ref) is \(o_p(\sqrt{n})\). This allows us to obtain doubly-robust inference for the CATE. \footnotetext{The requirement \(\lambda_{\alpha,j}/\lambda_{\gamma,j} \geq c\) may seem a bit unnatural, but it can be enforced in practice without upsetting any assumptions by setting the linear penalty \( \lambda_{\alpha,j}^{\text{\tiny ratio}} := \max\{\lambda_{\gamma,j}/5, \lambda_{\alpha, j}\} .\) In simulations, we find this constraint is rarely binding.}
Below, we provide some intuition for how this result is obtained and the role our particular estimating equations play in establishing this fact. We focus on control of the vector \(\vB^k\), defined in (ref), which measures the bias passed on from first-stage estimation to the second-stage estimate \(\widehat\beta^k\). Limiting the size of \(\vB^k\) is crucial in showing convergence of \(\widehat\beta^k\) to the true parameter \(\beta^k\) and thus consistency of the nonparametric estimator \(\widehat g(x)\).
For exposition, we consider a single term of (ref), \(\vB^k_j\), which roughly measures the first stage estimation bias taken on from adding the \(j^\text{th}\) basis term to our series approximation of \(g_0(x)\). The discussion that follows is a bit informal, instead of considering the derivatives with respect to the true parameters below our proof strategy will directly use the Kuhn-Tucker conditions of the optimization routines in (ref)-(ref). However, the general intuition is the same as is used in the proofs.
In addition to the doubly-robust identification property (ref), the aIPW signal is typically useful in the high-dimensional setting because it obeys an orthogonality condition at the true values \((\pi^\star,m^\star)\):\footnote{Robustness and orthogonality are indeed closely related, see Theorem 6.2 in nmf-1994 for a discussion.}
When both the propensity score model and outcome regression model are correctly specified we can (loosely speaking) examine the bias \(\vB_j^k\) by replacing \(\bar\pi_j = \pi^\star\) and \(\bar m_j = m^*\) and considering the following first order expansion:
By orthogonality of the aIPW signal the gradient term is close to zero, which guarantees that the bias is asymptotically negligible even if the nuisance parameters converge slowly to the true values, \(\pi^\star\) and \(m^\star\).\footnote{Typically all that is required is that \(\|\hat\pi_j - \pi^\star\| = o_p(n^{-1/4})\) and \(\|\hat m_j - m^\star\| = o_p(n^{-1/4})\) in order to make the second order remainder term \(\sqrt{n}\)-negligible} This allows the researcher to ignore first stage nuisance parameter estimation error and treat \(\pi^\star\) and \(m^\star\) as known when analyzing the asymptotic properties of the second stage series estimator. Indeed, since the aIPW signal orthogonality holds conditional on \(Z = (Z_1,X)\), if both models are correctly specified only a single pair of first stage estimators would be needed to provide control over all the elements in \(\vB^k\). This is the approach followed by SC-2020.
So long as either one of \(\bar\pi_j = \pi^\star\) or \(\bar m_j = m^\star\), double robustness of the aIPW signal (ref) still delivers identification: \(\E[p_j(X)Y_1] \approx \E_n[p_j(X)Y(\bar\pi_j,\bar m_j)\). However, the aIPW orthogonality tells us nothing about the expectation of the gradient away from the true parameters, \(\pi^\star, m^\star\); if either \(\bar\pi_j \neq \pi^\star\) or \(\bar m_j \neq m^\star\) there is no reason to believe that the gradient on the right hand side of (ref) is mean zero when evaluated instead at \(Y(\bar\pi_j,\bar m_j)\). In general, the bias \(\vB_j^k\) will then diminish at the rate of convergence of our nuisance parameters. Because we have high dimensional controls, this convergence rate will generally be much slower than the standard nonparametric rate Newey-1997,BCCK-2015.
To get around this, we design the first-stage objective functions (ref)-(ref) such that the resulting first-order conditions control the bias passed on to the second stage. Consider the following expansion instead around the limiting parameters \(\bar\gamma_j\) and \(\bar\alpha_k\).
After substituting the forms of \(\bar\pi_j(z) = \pi(z;\bar\gamma_j)\) and \(\bar m_j(z) = m(z;\bar\alpha_j)\) described in (ref) and differentiating with respect to \(\gamma_j\) and \(\alpha_j\) we obtain
However, by definition \(\bar\gamma_j\) and \(\bar\alpha_j\) solve the minimization problems defined in (ref)-(ref), the population analogs of our finite sample estimating equations. The first order conditions of these minimization problems yield
Examining the first order conditions in (ref), we see that they exactly give us control over the gradient (ref). Under suitable convergence of the first stage parameter estimates, this guarantees the bias examined in expansion (ref) is negligible even under misspecification of the propensity score or outcome regression models.
Control of this gradient under misspecification is not provided using other estimating equations, such as maximum likelihood for the logistic propensity score model or ordinary least squares for the linear outcome regression model. Moreover, control over the gradient of \(\vB_j^k\) from (ref) is not provided by the first-order conditions for \(\bar\gamma_l\) and \(\bar\alpha_l\) for \(l\neq j\):
Showing that the inference procedure of (ref) remains valid at all points \(x\in\calX\) under misspecification requires showing negligible first stage estimation bias for any linear combination of the vector (ref). As outlined above, this requires using \(k\) separate pairs of nuisance parameter estimator to obtain \(k\) separate pairs of first order conditions, one for each term of the vector.
In this section, we present the main consistency and distributional results for our second-stage estimator \(\widehat g(x)\) described in (ref). A full set of second stage results, including pointwise and uniform linearization lemmas and uniform convergence rates, can be found in (ref). The first set of results is established under the following condition, which limits the bias passed from first-stage estimation onto the second-stage estimator. In particular, (ref) implies that the bias vector \(\vB^k\) from (ref) satisfies \(\|\vB^k\| = o_p(n^{-1/2})\).
Via (ref) we can see that is a logistic propensity score model and a linear outcome regression model and estimating the first stage models using the estimating equations (ref)-(ref), (ref) can be achieved under (ref) and the sparsity bound
If the researcher were to assume different parametric forms for the first stage model, different first estimating equations would have to be used to obtain doubly-robust estimation and inference. However, so long as the (ref) can be established at the limiting values of the first stage models, the results of this section hold.
Having dealt with the first stage estimation error, the main complication remaining is that under misspecification the aIPW signals \(Y(\hat\pi_j, \hat m_j)\) for \(j = 1,\dots,k\) do not all converge to the same limiting values. However, so long as at least one of the first stage models is correctly specified, all of the limiting aIPW signals have the same conditional mean, \(g_0(x)\). In the standard setting, consistency of nonparametric estimator relies on certain conditions on the error terms. In our setting, we require that these assumptions hold uniformly over \(k\) the error terms. We note though that there is a non-trivial dependence structure between that limiting aIPW signals. This strong dependence gives plausibility to our uniform conditions. For example, if the logistic propensity score model is correctly specified and the limiting outcome regression models are uniformly bounded conditional on \(Z\), our conditions reduce exactly to the conditions of BCCK-2015. In general, however, the uniform conditions suggest that a degree of undersmoothing is optimal when implementing our estimation procedure.
Pointwise inference relies on the following assumption in tandem with (ref).
As mentioned, these are exactly the conditions required by BCCK-2015, with the modification that the bounds on conditional variance and other moment conditions on the error term hold uniformly over \(j = 1,\dots,k\). The assumptions on the series terms being used in the approximation can be shown to be satisfied by a number of commonly used functional bases, such as polynomial bases or splines, under adequate normalizations and smoothness of the underlying regression function. Readers should refer to Newey-1997, Chen-2007, or BCCK-2015 for a more in depth discussion of these assumptions.\footnote{In practice, we recommend the use of B-splines in order to to satisfy the first requirement that the basis functions are weakly positive and to reduce instability of the convex optimization programs described in (ref)-(ref).}
Under these assumptions, the variance of our second stage estimator is governed by one of the following variance matrices:
where \(\circ\) represents the Hadamard (element-wise) product and, abusing notation, for a vector \(a \in \SR^k\) and scalar \(c\in \SR\) we let \(a + c = (a_i + c)_{i=1}^k\). Later on, we establish the validity of the plug-in analog \(\hat\Omega\) (ref), as an estimator of these matrices.
(ref) shows that the estimator proposed in (ref) has a limiting gaussian distribution even under misspecification of either first stage model. This allows for doubly-robust pointwise inference after establishing a consistent variance estimator.
Next, we turn to strengthening the pointwise results to hold uniformly over all points \(x\in\calX\). This requires stronger conditions. we make the following assumptions on the tail behavior of the error terms which strengthens (ref).
As before, (ref) is very similar to its analogue in BCCK-2015, with the modification that the conditions are required to hold for \(\bar\eps_k\) as opposed to \(\eps_k\). Under this assumption, we derive doubly-robust uniform rates of convergence uniform inference procedures for the conditional counterfactual outcome \(g_0(x)\).
(ref) establishes conditions under which we obtain a doubly-robust strong approximation of the empirical process \(x \mapsto \sqrt{n}(\widehat g(x) - g_0(x))\) by a Gaussian process. After establishing consistent estimation of the matrix \(\Omega\), this strong approximation result allows us to show validity of the uniform confidence bands described in (ref). As noted by BCCK-2015, this is distinctly different from a Donsker type weak convergence result for the estimator \(\widehat g(x)\) as viewed as a random element of \(\ell^\infty(X)\). In particular, the covariance kernel is left completely unspecified and in general need not be well behaved.
We establish that the estimator \(\widehat\Omega\) proposed in (ref) is a consistent estimator of the true limiting variance \(\Omega\), where \(\Omega = \tilde\Omega\) in general but if \(\bar R_{2n} = o_p(a_n^{-1})\) then \(\Omega = \Omega_0\). To do so, we rely on the second stage assumptions (ref) as well as the following condition limiting the first stage estimation error passed on to the variance estimator \(\widehat\Omega\).
Via (ref) we can establish (ref) under (ref) as well as the additional sparsity bound\footnote{The sparsity bound (ref) required for consistent variance estimation can be significantly sharpened if the researcher is willing to use a cross fitting procedure, using one sample to estimate the nuisance parameters and another to evaluate the aIPW signal. This is because one could more directly follow SC-2020 and control alternate quantities with bounds that converge more quickly to zero.}
(ref) establishes that pointwise inference based on the test statistic described in (ref), obtained by replacing \(\Omega\) in (ref) with the consistent estimator \(\widehat\Omega\), is doubly-robust. Hypothesis tests based on the test statistic as well as pointwise confidence intervals for \(g_0(x)\) remain valid even if one of the first stage parameters is misspecified.
We now establish the validity of uniform inference based on the gaussian bootstrap critical values \(c_u^\star(1-\alpha)\) defined in (ref).
In conjunction with (ref), (ref) and (ref), (ref) shows the validity of the uniform inference procedure described in (ref).
Up to now, we have mainly focused on doubly-robust estimation and model-assisted inference for the function \[ g_0(x) = \E[Y_1 \mid X= x] .\] We conclude by noting that we can use a symmetric procedure to obtain model-assisted inference for the additional conditional counterfactual outcome \[ \tilde g_0(x) = \E[Y_0 \mid X = x] .\] To do so, we use the alternate aIPW signal \[ Y_0(\pi_0, m_0) = \frac{(1-D)Y}{1-\pi_0(Z)} + \left(\frac{1-D}{1-\pi_0(Z)} - 1\right)m_0(Z) \] where as before the true value for \(\pi^\star_0(z) = \Pr(D = 1\mid Z = z)\) but now \( m_0^\star(z) = \E[Y \mid D= 0, Z= z]\). To estimate these nuisance models we again assume a logistic form for the propensity score model \(\pi_0(z) = \pi(z; \gamma^0)\) and a linear form for the outcome regression model \(m_0(z) = m(z, \alpha^0)\) as in (ref) and use a separate estimation procedure for each basis term in our series approximation of \(\tilde g_0(x)\). The estimating equations we use to estimate each \(\gamma_j^0\) and \(\alpha_j^0\) differ from those in (ref)-(ref) however, and are instead given
which under the natural analog of (ref) converge uniformly to population minimizers:
Letting \(\bar\pi_{0,j}(z) = \pi(z, \bar\gamma^0_j)\), and \(\bar m_{0,j}(z) = m(z, \bar\alpha_j^0)\) we can repeat the decomposition of (ref), expressing \(\tilde Y(\bar\pi_{0,j}, \bar m_{0,j})\) as functions of the parameters \(\bar\gamma_j^0\) and \(\bar\alpha_j^0\) and show that the first order conditions for \(\bar\gamma_j^0\) and \(\bar\alpha_j^0\) directly control the bias passed on to the second stage nonparametric estimator for \(\tilde g_0(x)\). Convergence rates and validity of inference then follow from symmetric analysis of the results in (ref). Combining estimation and inference of the two conditional counterfactual outcomes then gives a doubly-robust estimator and inference procedure for the CATE. To perform inference on the CATE we can use the variance matrix \[ \bar\Omega = \Omega_0 + \Omega_1 - 2\Omega_2 \] where \(\Omega_0\) is as in (ref) but \(\Omega_1\) and \(\Omega_2\) are given
where \(\eps_{0,j}^k = Y_0(\bar\pi_{0,j},\bar m_{0,j}) - \tilde g_0(x)\) and \(\eps_0^k = (\eps_{0,1}^k,\dots,\eps_{0,k}^k)'\). These matrices can be consistently estimated using their natural empirical analogs as in (ref).
We investigate the finite-sample performance of the doubly-robust estimator and inference procedure via simulation study. We find that our proposed estimation procedure retains good coverage properties even under misspecification.
Observations are generated i.i.d. according to the following distributions The error term is generated following \(\epsilon \sim N(0, 1) \). The controls are set \(Z_i = (Z_{1i},X_i) \in \SR^{d_z}\) where \(d_z= 100\), \(X\sim U(1,2)\), and the independent regressors \(Z_1 \) are jointly centered Gaussian with a covariance matrix of the Toeplitz form
To capture misspecification, we let \(Z^\dagger\) be a transformation of the regressors in \(Z_1\) where \(Z_j^\dagger = Z_j + \max (0,1+Z_j)^2, \ \forall \ j=3,\dots,d_z\). Let sparsity control the number of regressors in \(Z = (Z_1,X)\) entering the DGP.
where the constants \(p_1\) and \(p_2\) differ in various simulation setups but are always set so that the average probability of treatment is about one half. To consider various degrees of high-dimensionality, we implement \( N \in \{500, 1000\} \) with \(d_z = 100\). For (S1), sparsity\(=6\); for (S2), sparsity\(=4\); and, for (S3), sparsity\(=5\). Results are reported for \(S=1,000\) repeated simulations.
To select the first stage penalty parameters, we implement the multiplier bootstrap procedure described in (ref). The constants \(c_{\gamma,j}\) and \(c_{\alpha,j}\) in the pilot penalty parameters (ref) are selected via cross validation from a set of size 5. To select the final bootstrap penalty parameter we set \(c_0 = 1.1\) and select the \(95^\text{\tiny th}\) quantile of \(B=10000\) bootstrap replications. In our second-stage estimation, we use a b-spline basis of size \(k=3\). B-splines are implemented from the R package splines2 splines2-paper, which uses the specification detailed in perperoglou2019review. In the tables below, we refer to our method as MA-DML (model assisted double machine learning).
We compare our proposed estimator and inference procedure to that of SC-2020, which projects a single aIPW signal onto a growing series of basis terms. In implementing this DML method, we use the standard \(\ell_1\)-penalized maximum likelihood (MLE) and ordinary least squares (OLS) loss functions to estimate the first stage propensity score and outcome regression models, respectively.\footnote{Vira Semenova provides several example R scripts implementing DML: \url{https://sites.google.com/view/semenovavira/research}.}
Estimation error is studied for the target parameter \(g_0 (x)= \E [ Y| D=1, X=x]\) over a grid of 100 points spaced across \(x\in[1,2]\), i.e. the support of \( X \). We study average coverage across simulations of each method's pointwise (at \(x = 1.5 \)) and uniform confidence intervals. To compare the estimation error for the target parameter \( g(x) \) across the two different estimators \( \widehat g_s (x) \) for each simulation \( s = 1,\dots, S \), we utilize integrated bias, variance, and mean-squared error where \( \Bar{g} (x) = S^{-1} \sum^S_{s=1} \widehat g_s (x), \)
Table (ref) presents the simulation results for all three specifications (S1)-(S3) for \(n=500\) and \(n = 1000\). Integrated squared bias, variance, and mean squared error are presented in columns (1)-(3), respectively. Pointwise and uniform coverage results are presented in columns (4)-(7).
For pointwise and uniform coverage under correct specification regime (S1), MA-DML has some slight improvements. Under misspecification DGPs (S2) and (S3), the pointwise coverage of MA-DML is closer to the targets except in the $N=1000$ and (S2) case where it slightly underperforms. However, MA-DML has a notable improvement over DML in the (S3) case when $N=1000.$ Similarly, MA-DML outperforms DML in three of the four misspecified regimes, i.e. all but (S3) when $N=500$ where \textit{MA-DML} has over-coverage. Under (S2) when $N=1000,$ both methods are markedly deterioated uniform coverage, although \textit{MA-DML} is noticably closer to target.
In regards to estimation error, in four of the six settings, MA-DML has a lower MSE than DML where regardless of sample size MA-DML underperforms in (S3). Notably, it does appear MA-DML has substantially smaller IBias$^2$ across the DGPs.
Finally, we were surprised to find for both estimators that coverage properties, in general, improve under the higher-dimensional regime of $N=500$ with $d_z=100$ compared to $N=1,000$ and $d_z=100.$ In particular, with a higher ratio of covariates to observations, the uniform coverage properties under regime (S2) were substantially better. The estimation error results were in line with our priors as the higher-dimensional regime sees in general higher estimation errors for both methods.
For coverage under correct specification, we did anticipate the underperformance of MA-DML given it is designed to handle misspecification with the cost of other estimators outperforming under correct specification. Additionally, we attribute the poor uniform coverage in DGP (S2) for both estimators under $N=1,000$ to a lack of a rich enough cross-validation given the performance was improved under a more difficult regime when the number of observations drops to $N=500.$ The integrated bias of MA-DML is lower across the various DGPs compared to DML. Following the discussion in (ref) this is expected since the first stage estimating equations for the model assisted procedure are specifically designed to minimize the bias passed on to the second stage estimator. However, the model assisted procedure has higher values of integrated variance compared to the standard procedure, which could be attributable to the use of \(k\) distinct first-stage estimations.
Our findings should not be interpreted as a critique of the SC-2020 benchmark method, whose work we rely on and were inspired by.
We apply the model assisted estimator to estimate the effect of maternal smoking on infant birthweight conditional on the age of the mother. We use the Cattaneo_2010 dataset which can be found online on the Stata website.\footnote{The dataset can be downloaded \href{http://www.stata-press.com/data/r13/cattaneo2.dta}{here}.} The dataset describes each infant's birthweight in grams, \(Y\), whether or not the mother smoked during pregnancy, \(D=1\) indicating smoking, and a number of covariates containing information on the mother's health and socioeconomic background, \(Z = (X,Z_1)\), where \(X\) represents the conditioning variable, maternal age. A full summary of the data used as well as additional details/analysis from our empirical analysis can be found in (ref).
We compare the model assisted estimator of the CATE against one where standard MLE and OLS loss functions are used to estimate the first stage propensity score and outcome regression models. We also qualitatively compare our results to Zimmert_Lechner_CATE, who use a kernel based approach to estimate the CATE in this setting. While this sort of comparison is not perfect since we do not know the true DGP, this setting is advantageous for analysis since we strongly expect that (i) the effect of smoking on birthweight will be negative and (ii) this effect should grow stronger in magnitude as the age of the mother increases. These hypotheses have been corroborated by other work that examines the conditional average treatment effect in this setting Zimmert_Lechner_CATE,Abreya_2006,Lee_Ryo_Wang_2016.
(ref) displays our main results from implementing both the model assisted and standard MLE/OLS estimation procedures. After removing the top 3% and bottom 3% of smoker and non-smoker birthweights by maternal age, we select the penalty parameters for the first stage models via the bootstrap procedure described in (ref). The pilot penalty parameters are uniformly taken to be equal to zero, so that the residuals used in the bootstrap procedure are generated from non-regularized estimations. We take \(c_0 = 2\) in (ref) and and select the first stage penalty parameters using the 99\textsuperscript{th}, 95\textsuperscript{th}, and 90\textsuperscript{th} quantiles of the bootstrap distribution. For the second stage basis functions we implement second degree b-splines with 3 knots via the splines2 package in R splines2-paper.
Consistent with prior work, both estimators of the CATE suggest that the effect of smoking on birthweight becomes more negative with age. Both estimation procedures also generally produces negative estimates for the CATE, but it should be noted that for the lowest levels of penalization the model assisted CATE estimate suggests a slightly positive effect of smoking for particularly young mothers, though this difference is not significantly different from zero. The shapes of the estimated functions remain relatively stable under various sizes of the penalty parameter, though the model assisted procedure displays a bit more sensitivity to the level of regularization introduced.\footnote{Numerically solving the minimization problems in (ref)-(ref) also typically requires more iterations to converge than solving the standard MLE/OLS minimization problems.}
For the most part, the effects found here are similar to those found in Zimmert_Lechner_CATE, though the effects estimated using standard first stage loss functions have somewhat larger magnitudes and in general both series estimation procedures seem to give less reasonable results on the boundaries. An advantage of using a series second stage however, compared to the kernel first stage of Zimmert_Lechner_CATE, is the existence of the uniform confidence bands displayed. Reassuringly, the estimates of Zimmert_Lechner_CATE seem to be within the 95% uniform confidence bands generated by the model assisted estimator.
As a robustness check, we also try estimating the treatment effect using first degree b-splines instead of second degree splines. These results are displayed in (ref). Again, we find that the effect of smoking on child birthweight is almost uniformly negative regardless of estimation procedure used or choice of penalty parameter. The shape of the estimated CATE function using a standard MLE/OLS first stage is very stable to penalty choice here while the shape of the model assisted CATE function displays a bit more instability here at the two lower levels of regularization.
Finally, (ref) reports the smoothed average treatment effect estimates taken from averaging the model assisted CATE estimates from (ref) across observations. Again, these estimates are generally in line with prior work
Estimation of conditional average treatment effects with high dimensional controls typically relies on first estimating two nuisance parameters: a propensity score model and an outcome regression model. In a high-dimensional setting, consistency of the nuisance parameter estimators typically relies on correctly specifying their functional forms. While the resulting second-stage estimator for the conditional average treatment effect typically remains consistent even if one of the nuisance parameters is inconsistent, the confidence intervals may no longer be valid.
In this paper, we consider estimation and valid inference on the conditional average treatment effect in the presence of high dimensional controls and nuisance parameter misspecification. We present a nonparametric estimator for the CATE that remains consistent at the nonparametric rate, under slightly modified conditions, even under misspecification of either the logistic propensity score model or linear outcome regression model. The resulting Wald-type confidence intervals based on this estimator also provide valid asymptotic coverage under nuisance parameter misspecification.
\singlespacing