The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
77,223 characters
Doubly-Robust Inference for Conditional Average Treatment Effects with High-Dimensional Controls
\begin{doublespace}
\maketitle
\end{doublespace}
\begin{abstract}
Plausible identification of conditional average treatment effects (CATEs) may rely on controlling for a large number of variables to account for confounding factors. In these high-dimensional settings, estimation of the CATE requires estimating first-stage models whose consistency relies on correctly specifying their parametric forms. While doubly-robust estimators of the CATE exist, inference procedures based on the second stage CATE estimator are not doubly-robust. Using the popular augmented inverse propensity weighting signal, we propose an estimator for the CATE whose resulting Wald-type confidence intervals are doubly-robust. We assume a logistic model for the propensity score and a linear model for the outcome regression, and estimate the parameters of these models using an \(\ell_1\) (Lasso) penalty to address the high dimensional covariates. Our proposed estimator remains consistent at the nonparametric rate and our proposed pointwise and uniform confidence intervals remain asymptotically valid even if one of the logistic propensity score or linear outcome regression models are misspecified.
\end{abstract}
\newpage
\section{Introduction}
\label{sec:introduction}
Consider a potential outcomes framework \citep{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):
\begin{equation}
\label{eq:CATE-def}
\E[Y_1 - Y_0 \mid X = x]
.\end{equation}
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 \citet{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 \citep{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 \citet{Tan-2018} wherein doubly-robust inference is developed for the average treatment effect. Following \citet{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 \citet{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 \citep{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 \citet{BCCK-2015}, though in certain special cases collapse exactly to the conditions of \citet{BCCK-2015}.
\paragraph{Prior Literature.}
\citet{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. \citet{SC-2020} take advantage of new results for series estimation in \citet{BCCK-2015} and consider series estimation of functional target parameters after high-dimensional nuisance estimation.\footnote{\citet{fan2022estimation} provides a similar analysis using a second stage kernel estimator.}
In the same setting as this paper, \citet{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, \citet{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. \citet{SRR-2019} extends the analysis of \citet{Tan-2018} to consider doubly-robust inference for a larger class of finite dimensional target parameters with bilinear influence functions.
\citet{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.
\citet{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 \citep{tibshirani1996regression,BRT-2009,BC-2013,CLC-2021-CVLasso} and \(\ell_1\)-regularized M-estimation in high dimensional settings \citep{vanDerGreer2016,Tan-2017}.
\paragraph{Paper Structure.} This paper proceeds as follows. Section~\ref{sec:setup} defines the problem and introduces our methods for estimation and inference. Section~\ref{sec:theory-overview} 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{sec:first-stage} 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{sec:cate-wrapup} ties up a technical detail. Section~\ref{sec:simulations} provides evidence from a simulation study while Section~\ref{sec:empirical} applies our proposed estimator to examine the effect of maternal smoking on infant birth weight. Section~\ref{sec:conclusion} 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\).
\section{Setup}
\label{sec:setup}
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:
\begin{equation}
\label{eq:target-paramter}
g_0(x) := \E[Y_1\mid X= x]
.\end{equation}
Doubly-robust estimation and inference on the other conditional counterfactual outcome, \(\E[Y_0\,|\,X=x]\), follows a similar procedure and is described in \Cref{sec:cate-wrapup}. The procedures can be combined for doubly-robust estimation and inference for the CATE.
\subsection{Setting}
\label{subsec:setting}
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)\).
\begin{assumption}[Identification]
\label{assm:identification}
\leavevmode
\begin{enumerate}[(i)]
\item \(\{Y_i,D_i,Z_i\}_{i=1}^n\) are independent and identically distributed.
\item \((Y_1,Y_0)\perp D \mid Z\).
\item There exists a value \(\eta \in (0,1) \) such that \(\eta < \E[D\mid Z= z] < 1-\eta\) almost surely in \(Z\).
\end{enumerate}
\end{assumption}
To obtain doubly-robust estimation and inference we use the augmented inverse propensity weighted (aIPW) signal,
\begin{equation}
\label{eq:aIPW-signal}
Y(\pi, m) = \frac{DY}{\pi(Z)} - \left(\frac{D}{\pi(Z)} - 1\right)m(Z),
\end{equation}
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 \Cref{assm:identification}, 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\),
\begin{equation}
\label{eq:double-robustness}
\begin{split}
\E[Y_1 \mid X = x] &= \E[Y(\pi^\star, m^\star) \mid X = x] \\
&= \E[\,Y(\pi, m^\star)\;\mid X = x] \\
&= \E[\,Y(\pi^\star, m)\; \mid X = x].
\end{split}
\end{equation}
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:
\begin{equation}
\label{eq:basis-terms}
p^k(x) := \left(p_1(x),\dots,p_k(x)\right)' \in \SR_+^k
.\end{equation}
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\):
\begin{equation}
\label{eq:second-stage-setup}
\begin{split}
Y(\bar\pi_j, \bar m_j)
&= g_0(x) + \eps_j,\,\;\;\;\;\;\;\;\;\;\;\;\;\;\E[\eps_j\mid X] = 0\\
&= g_k(x) + r_k(x) + \eps_j
\end{split}
\end{equation}
where \(g_0(x)\) is the conditional counterfactual outcome~\eqref{eq:target-paramter}, \(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 \eqref{eq:second-stage-setup}, 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)\):
\begin{equation}
\label{eq:beta-k-population}
\begin{split}
\beta^k &:= Q^{-1}\E[p^k(X)Y_1] \\
&\;= Q^{-1}\E[p^k(X)Y(\pi^\star, m^\star)] \\
&\;= Q^{-1}\E[p^k(X)Y(\bar\pi_j,\bar m_j)],\;\; \forall j = 1,\dots,k.
\end{split}
\end{equation}
\footnotetext{\Cref{sec:first-stage-second-stage} 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.}
\subsection{Estimator and Inference Procedure}
\label{subsec:estimator-inference}
We assume a logistic regression form for the propensity score model and a linear form for the outcome regression model:
\begin{equation}
\label{eq:nuisance-parameter-functional-forms}
\begin{split}
\pi(Z; \gamma) &= \left(1 + \exp(-\gamma'Z)\right)^{-1}, \\
m(Z; \alpha) &= \alpha'Z.
\end{split}
\end{equation}
For each \(j=1,\dots,k,\) the parameters of \eqref{eq:nuisance-parameter-functional-forms}, \(\gamma,\alpha \in \SR^{d_z} ,\) are estimated by
\begin{align}
\label{eq:gamma-j-estimating-equation}
\widehat\gamma_j &:= \arg\min_\gamma\, \E_n[p_j(X)\{De^{-\gamma'Z} + (1-D)\gamma'Z\}] + \lambda_{\gamma,j}\|\gamma\|_1, \\
\label{eq:alpha-j-estimating-equation}
\widehat\alpha_j &:= \arg\min_\alpha\, \E_n[p_j(X)De^{-\widehat\gamma_j'Z}(Y - \alpha'Z)^2]/2 + \lambda_{\alpha,j}\|\alpha\|_1.
\end{align}
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 \Cref{sec:theory-overview}. Under standard assumptions the parameter estimators \(\widehat\gamma_j, \widehat\alpha_j\) will converge uniformly over \(j=1,\dots,k\) to population minimizers
\begin{align}
\label{eq:gamma-bar-j}
\bar\gamma_j &:= \arg\min_\gamma\, \E[p_j(X)\{De^{-\gamma'Z} + (1-D)\gamma'Z\}], \\
\label{eq:alpha-bar-j}
\bar\alpha_j &:= \arg\min_\alpha \E[p_j(Z)De^{-\bar\gamma_j'Z}(Y - \alpha'Z)^2].
\end{align}
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
\begin{equation}
\label{eq:beta-k}
\widehat\beta^k = \widehat{Q}^{-1}\E_n
\begin{bmatrix} p_1(X)Y(\widehat\pi_1, \widehat m_1) \\ \vdots \\ p_k(X)Y(\widehat\pi_k, \widehat m_k) \end{bmatrix},
\end{equation}
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
\begin{equation}
\label{eq:omega-hat-definitions}
\widehat\Omega := \widehat Q^{-1}\E_n[\{p^k(X)\circ\widehat\eps^k\}\{p^k(X)\circ\widehat\eps^k\}' ]\widehat Q^{-1},
\end{equation}
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
\begin{equation}
\label{eq:confidence-bands}
\left[\underline i(x), \bar i(x)\right] := \left[\widehat g(x) - c^\star\left(1-\eta/2\right)\widehat\sigma(x),\; \widehat g(x) + c^\star\left(1-\eta/2\right)\widehat\sigma(x)\right]
.\end{equation}
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)\). \Cref{sec:theory-overview,sec:first-stage} show that, under standard sparsity and moment conditions, these pointwise and uniform inference procedures remain valid even under misspecification of either first-stage model.
\subsection{Penalty Parameter Selection}
\label{subsec:additional}
To select the penalty parameters \(\lambda_{\gamma,j}\) and \(\lambda_{\alpha,j}\) in \eqref{eq:gamma-j-estimating-equation}-\eqref{eq:alpha-j-estimating-equation} we propose a data driven two-step procedure based on the work of \citet{CS-2021}.
For each \(j= 0,1\dots,k,\) we start with pilot penalty parameters given by
\begin{equation}
\label{eq:pilot-penalty}
\lambda^{\text{\tiny pilot}}_{\gamma, j} = c_{\gamma,j}\times \sqrt{\frac{\ln^3(d_z)}{n}} \andbox \lambda_{\alpha,j}^{\text{\tiny pilot}} = c_{\alpha,j}\times \sqrt{\frac{\ln^3(d_z)}{n} }
\end{equation}
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 \Cref{assm:logistic-model-convergence}(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 \eqref{eq:gamma-j-estimating-equation}-\eqref{eq:alpha-j-estimating-equation} 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
\begin{equation}
\label{eq:residual-estimates}
\begin{split}
\widehat U_{\gamma,j} &:= -p_j(X)\{De^{-\widehat\gamma_j^{\text{\tiny pilot}'}Z} + (1-D)\} \\
\widehat U_{\alpha,j} &:= p_j(X)De^{-\widehat\gamma_j^{\text{\tiny pilot}'}Z}(Y - \widehat\alpha_j^{\text{\tiny pilot}'}Z).
\end{split}
\end{equation}
We then use a multiplier bootstrap procedure to select our final penalty parameters \(\lambda_{\gamma, j}\) and \(\lambda_{\alpha, j}\).
\begin{equation}
\label{eq:final-penalty-parameters}
\begin{split}
\lambda_{\gamma, j} &= c_0\times (1-\eps)\text{-quantile of} \max_{1\leq l \leq d_z} |\E_n[e_i\widehat U_{\gamma,j}Z_l]|\text{ given }\{Y_i, D_i, Z_i\}_{i=1}^n, \\
\lambda_{\alpha, j} &= c_0\times(1-\eps)\text{-quantile of} \max_{1\leq l \leq d_z} |\E_n[e_i\widehat U_{\alpha,j}Z_l]|\text{ given }\{Y_i, D_i, Z_i\}_{i=1}^n
\end{split}
\end{equation}
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 \Cref{thm:first-stage-convergence}. 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 \eqref{eq:final-penalty-parameters} 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 \eqref{eq:pilot-penalty} 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. \Cref{sec:alt-penalty-CV} provides details for this implementation as well as formally shows the modified conditions needed.
\section{Theory Overview}
\label{sec:theory-overview}
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 \citep{BCH-2013,Tan-2018} and is achieved under standard conditions in the \(\ell_1\)-regularized estimation literature \citep{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 \Cref{subsec:first-stage-intuition} 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.
\subsection{Uniform First-Stage Convergence}
\label{subsec:uniform-first-stage-convergence}
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:
\begin{assumption}[First Stage Convergence]
\label{assm:logistic-model-convergence}
\leavevmode
\begin{enumerate}[(i)]
\item The regressors \(Z\) are bounded, \(\max_{1 \leq l \leq d_z} |Z_l| \leq C_0\) almost surely.
\item The errors \(Y_1 - \bar m_j(Z)\) are uniformly subgaussian conditional on \(Z\) in the following sense. There exist fixed positive constants \(G_0\) and \(G_1\) such that for any \(j\):
\[
G_0 \E\left[\exp\big(\{Y_1 - \bar m_j(Z)\}^2/G_0^2\big) - 1\mid Z\right] \leq G_1^2
\]
almost surely.
\item There is a constant \(B_0\) such that \(\bar\gamma_j'Z \geq B_0\) almost surely for all \(j\).
\item There exist fixed constants \(\xi_0 > 1\) and \(1> \nu_0 > 0\) such that for each \(j = 1,\dots,k\) the following empirical compatability condition holds for the empirical hessian matrix \(\tilde\Sigma_{\gamma, j} := \E_n[De^{-\bar\gamma_j'Z}ZZ']\). For any \(b\in\SR^{d_z}\) and \(\calS_j = \{l: \bar\gamma_l\vee\bar\alpha_l\neq 0\}\):
\[
\sum_{l\not\in\calS_j} |b_l|\leq \xi_0\sum_{l\in\calS_j} |b_l|
\implies \nu_0^2\Big(\sum_{l\in\calS_j}|b_l|\Big)^2 \leq |\calS| \left(b'\tilde\Sigma_{\gamma, j}b\right)
.\]
\item There exist fixed constants \(c_u\) and \(C_U > 0\) such that for all \(j \leq k\), \(\E[U_{\gamma,j}^4] \leq (\xi_{k,\infty}C_U)^4\) and \(\min_{1\leq l\leq d_z}\E[U_{\gamma, j}^2 Z_l^2] \geq c_u\).
\item The constant \(\underline c_n\) is chosen such that \(\xi_{k.\infty} \lesssim \underline c_n\) and the following sparsity bounds hold for \(s_k = \max_{1 \leq j \leq k}|\calS_j|\)
\[
\frac{\xi_{k,\infty}s_k^2\bar c_n^2\ln^5(d_zn)}{n} \to 0,\andbox \frac{\xi_{k,\infty}^4\ln^7(d_zkn)}{n}\to 0
.\]
\end{enumerate}
\end{assumption}
The first part of \Cref{assm:logistic-model-convergence} 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 \citep{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 \Cref{assm:logistic-model-convergence} 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 \citet{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 \eqref{eq:gamma-j-estimating-equation}-\eqref{eq:alpha-j-estimating-equation} as well as some assumptions on the order of the constants \(c_{\gamma,j}\) and \(c_{\alpha,j}\) in \eqref{eq:pilot-penalty}.
\begin{lemma}[First-Stage Convergence]
\label{thm:first-stage-convergence}
Suppose that \Cref{assm:logistic-model-convergence} holds. In addition assume that \(c_0 > (\xi_0 + 1)/(\xi_0-1)\), \(k/n \to 0\), \(k\eps \to 0\), and there is a fixed constant \(c > 0\) such that for all \(j\), \(\lambda_{\alpha, j}/\lambda_{\gamma,j} \geq c\).\footnotemark Then the following weighted means converge uniformly in absolute value at least at rate:
\begin{equation}
\label{eq:uniform-mean-convergence}
\max_{1\leq j\leq k}\left|\E_n[p_j(X)Y(\widehat\pi_j, \widehat m_j)] - \E_n[p_j(X)Y(\bar\pi_j, \bar m_j)]\right|
\lesssim_P
\frac{s_k\,\xi_{k,\infty}^2\ln(d_z)}{n}
\end{equation}
and in empirical mean square at least at rate:
\begin{equation}
\label{eq:second-moment-convergence}
\max_{1\leq j\leq k}\E_n[p_j^2(X)(Y(\widehat\pi_j,\widehat m_j) - Y(\bar\pi_j,\bar m_j))^2] \lesssim_P \frac{s_k^2\,\xi_{k,\infty}^4 \ln(d_z)}{n}
\end{equation}
\end{lemma}
\Cref{thm:first-stage-convergence} 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 \eqref{eq:uniform-mean-convergence} and \eqref{eq:second-moment-convergence} 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.}
\subsection{Managing First-Stage Bias}
\label{subsec:first-stage-intuition}
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 \eqref{eq:beta-means}, 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)\).
\begin{equation}
\label{eq:beta-means}
\vB^k :=\E_n\begin{bmatrix} p_1(X)\left\{Y(\widehat\pi_1, \widehat m_1) - Y(\bar\pi_1, \bar m_1)\right\} \\ \vdots \\ p_k(X)\left\{Y(\widehat \pi_k, \widehat m_k) - Y(\bar \pi_k, \bar m_k)\right\}\end{bmatrix}.
\end{equation}
For exposition, we consider a single term of \eqref{eq:beta-means}, \(\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 \eqref{eq:gamma-j-estimating-equation}-\eqref{eq:alpha-j-estimating-equation}. However, the general intuition is the same as is used in the proofs.
In addition to the doubly-robust identification property \eqref{eq:double-robustness}, 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 \citet{nmf-1994} for a discussion.}
\begin{equation}
\label{eq:neyman-orthogonality}
\E[\nabla_{\pi,m} Y(\pi^\star, m^\star) \mid Z] = 0
.\end{equation}
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:
\begin{equation}
\label{eq:taylor-expansion-true-parameters}
\begin{split}
\vB_j^k=
\E_n[p_j(X)Y(\widehat \pi_j, \widehat m_j)]
&- \E_n[p_j(X)Y(\pi^\star, m^\star)] \\
&= \underbrace{\E_n[p_j(X)\nabla_{\pi, m}\,Y(\pi^\star, m^\star)]}_{O_p(n^{-1/2})\text{ by \eqref{eq:neyman-orthogonality}}}\begin{bmatrix} \widehat \pi_j - \pi^\star \\ \widehat m_j - m^\star \end{bmatrix} + o_p(n^{-1/2}).
\end{split}
\end{equation}
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 \citet{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~\eqref{eq:double-robustness} 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 \eqref{eq:taylor-expansion-true-parameters} 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 \citep{Newey-1997,BCCK-2015}.
To get around this, we design the first-stage objective functions \eqref{eq:gamma-j-estimating-equation}-\eqref{eq:alpha-j-estimating-equation} 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\).
\begin{equation}
\label{eq:taylor-expansion-gamma-alpha}
\begin{split}
\vB_j^k = \E_n[p_j(X)Y(\widehat \pi_j, \widehat m_j)]
&- \E_n[p_j(X)Y(\bar\pi_j, \bar m_j)] \\
&= \E_n[p_j(X)\nabla_{\gamma_j, \alpha_j}\,Y(\bar\pi_j, \bar m_j)]
\begin{bmatrix} \widehat \gamma_j - \bar\gamma_j \\ \widehat \alpha_j - \bar\alpha_j \end{bmatrix} + o_p(n^{-1/2})
\end{split}
\end{equation}
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 \eqref{eq:nuisance-parameter-functional-forms} and differentiating with respect to \(\gamma_j\) and \(\alpha_j\) we obtain
\begin{equation}
\label{eq:gradient-plim-parameters}
\E[p_j(X)\nabla_{\gamma_j,\alpha_j}\,Y(\bar\pi_j,\bar m_j)]
= \E\begin{bmatrix}
-p_j(X)De^{-\bar\gamma'Z}(Y- \bar\alpha'Z)Z \\
p_j(x)\{D(1+e^{-\bar\gamma'Z})Z + Z\}
\end{bmatrix}
\end{equation}
However, by definition \(\bar\gamma_j\) and \(\bar\alpha_j\) solve the minimization problems defined in \eqref{eq:gamma-bar-j}-\eqref{eq:alpha-bar-j}, the population analogs of our finite sample estimating equations. The first order conditions of these minimization problems yield
\begin{equation}
\label{eq:foc-plim-parameters}
\begin{split}
\E
\underbrace{\overbrace{\begin{bmatrix}
p_j(X)\{D(1 + e^{\bar\gamma'Z})Z + Z\} \vspace{0.2cm}\\
p_j(X)De^{-\bar\gamma'Z}(DY - \bar\alpha'Z)Z
\end{bmatrix}}^{\text{First order condition of \(\bar\gamma_j\)}}}_{\text{First order condition of \(\bar\alpha_j\)}}
= 0 \;\implies\; \E[p_j(X)\nabla_{\gamma_j,\alpha_j}Y(\bar\pi_j,\bar m_j)] = 0
\end{split}
\end{equation}
Examining the first order conditions in \eqref{eq:foc-plim-parameters}, we see that they exactly give us control over the gradient \eqref{eq:gradient-plim-parameters}. Under suitable convergence of the first stage parameter estimates, this guarantees the bias examined in expansion \eqref{eq:taylor-expansion-gamma-alpha} 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 \eqref{eq:beta-means} is not provided by the first-order conditions for \(\bar\gamma_l\) and \(\bar\alpha_l\) for \(l\neq j\):
\begin{equation}
\label{eq:foc-l-does-not-control-gradient-j}
\begin{split}
\E[p_j(X)\nabla_{\gamma_j, \alpha_j}Y(\bar\pi_j, \bar m_j)]
&= \E\begin{bmatrix}
-p_j(X)De^{-\bar\gamma'Z}(Y - \bar\alpha'Z)Z \vspace{0.2cm} \\
p_j(X)\{D(1 + e^{\bar\gamma'Z})Z + Z\}
\end{bmatrix} \\
&\neq \E
\underbrace{\overbrace{\begin{bmatrix}
p_l(X)\{D(1 + e^{\bar\gamma'Z})Z + Z\} \vspace{0.2cm}\\
p_l(X)De^{-\bar\gamma'Z}(Y - \bar\alpha'Z)Z
\end{bmatrix}}^{\text{First order condition of \(\bar\gamma_l\)}}}_{\text{First order condition of \(\bar\alpha_l\)}}
.\end{split}
\end{equation}
Showing that the inference procedure of \Cref{sec:setup} remains valid at all points \(x\in\calX\) under misspecification requires showing negligible first stage estimation bias for any linear combination of the vector~\eqref{eq:beta-means}. 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.
\section{Main Results}
\label{sec:first-stage}
In this section, we present the main consistency and distributional results for our second-stage estimator \(\widehat g(x)\) described in \Cref{sec:setup}. A full set of second stage results, including pointwise and uniform linearization lemmas and uniform convergence rates, can be found in \Cref{sec:additional-second-stage}. 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, \Cref{cond:no-effect} implies that the bias vector \(\vB^k\) from \eqref{eq:beta-means} satisfies \(\|\vB^k\| = o_p(n^{-1/2})\).
\begin{condition}[No Effect of First-Stage Bias]
\label{cond:no-effect}
\begin{equation}
\label{eq:no-effect-first-stage-bias}
\max_{1 \leq j \leq k}\big|\E_n[p_j(X)Y(\widehat\pi_j, \widehat m_j)] - \E_n[p_j(X)Y(\bar\pi_j, \bar m_j)]\big| = o_p(n^{-1/2}k^{-1/2})
.\end{equation}
\end{condition}
Via \Cref{thm:first-stage-convergence} 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 \eqref{eq:gamma-j-estimating-equation}-\eqref{eq:alpha-j-estimating-equation}, \Cref{cond:no-effect} can be achieved under \Cref{assm:logistic-model-convergence} and the sparsity bound
\begin{equation}
\label{eq:first-stage-sparsity-bound}
\frac{s_k\,\xi_{k,\infty}^2k^{1/2}\ln(d_z)}{\sqrt{n}} \to 0
.\end{equation}
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 \Cref{cond:no-effect} 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 \citet{BCCK-2015}. In general, however, the uniform conditions suggest that a degree of undersmoothing is optimal when implementing our estimation procedure.
\subsection{Pointwise Inference}
\label{subsec:pointwise-inference}
Pointwise inference relies on the following assumption in tandem with \Cref{cond:no-effect}.
\begin{assumption}[Second-Stage Pointwise Assumption]
\label{assm:second-stage-assumptions}
Let \(\bar \eps_k := \max_{1\leq j\leq k}|\eps_j|\). Assume that
\begin{enumerate}[(i)]
\item Uniformly over all \(n\), the eigenvalues of \(Q = \E[p^k(x)p^k(x)']\) are bounded from above and away from zero.
\item The conditional variance of the error terms is uniformly bounded in the following sense. There exist constants \(\underline\sigma^2\) and \(\bar\sigma^2\) such that for any \(j=1,2\dots\) we have that \(\underline{\sigma}^2 \leq \Var(\eps_j \mid X)\leq \bar\sigma^2 < \infty;\)
\item For each \(n\) and \(k\) there are finite constants \(c_k\) and \(\ell_k\) such that for each \(f\in \calG\)
\[
\|r_k\|_{L,2} = (\E[r_k(x)^2])^{1/2}\leq c_k\andbox \|r_k\|_{L,\infty} = \sup_{x\in\calX}|r_k(x)| \leq \ell_kc_k
.\]
\item \(\sup_{x\in\calX}\E[\bar\eps_k^2\,\bm{1}\{\bar\eps_k + \ell_kc_k > \delta\sqrt{n}/\xi_k\} \mid X= x]\to 0\) as \(n\to\infty\) and \(\sup_{x\in\calX}\E[\ell_k^2c_k^2\bm{1}\{\bar\eps_k + \ell_kc_k > \delta\sqrt{n}/\xi_k\} \mid X= x]\to 0\) as \(n\to\infty\) for any \(\delta > 0\).
\end{enumerate}
\end{assumption}
As mentioned, these are exactly the conditions required by \citet{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 \citet{Newey-1997}, \citet{Chen-2007}, or \citet{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 \eqref{eq:gamma-j-estimating-equation}-\eqref{eq:alpha-j-estimating-equation}.}
Under these assumptions, the variance of our second stage estimator is governed by one of the following variance matrices:
\begin{equation}
\label{eq:omega-definitions}
\begin{split}
\tilde\Omega &:= Q^{-1}\E[\{p^k(x)\circ(\eps^k + r_k)\}\{p^k(x)\circ(\eps^k + r_k)\}' ]Q^{-1} \\
\Omega_0 &:= Q^{-1}\E[\{p^k(x)\circ\eps^k\}\{p^k(x)\circ\eps^k\}']Q^{-1}\\
\end{split}
\end{equation}
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\) \eqref{eq:omega-hat-definitions}, as an estimator of these matrices.
\begin{theorem}[Pointwise Normality]
\label{thm:pointwise-normality}
Suppose that \Cref{cond:no-effect} and \Cref{assm:second-stage-assumptions} hold. In addition suppose that \(\xi_k^2\log k/n\to 0\). Then so long as either the logistic propensity score model or linear outcome regression model is correctly specified, for any \(\alpha\in S^{k-1}\):
\begin{equation}
\label{eq:thm-pointwise-normality-1}
\sqrt{n}\frac{\alpha'(\widehat\beta^k-\beta^k)}{\|\alpha'\Omega^{1/2}\|} \to_d N(0,1)
\end{equation}
where generally \(\Omega = \tilde\Omega\) but if \(\ell_kc_k \to 0\) then we can set \(\Omega = \Omega_0\). Moreover, for any \(x\in\calX\) and \(s(x) := \Omega^{1/2}p^k(x)\),
\begin{equation}
\label{eq:thm-pointwise-normality-2}
\sqrt{n}\frac{p^k(x)'(\widehat\beta^k - \beta^k)}{\|s(x)\|} \to_d N(0,1)
\end{equation}
and if the approximation error is negligible relative to the estimation error, namely \(\sqrt{n}r_k(x) = o(\|s(x)\|)\), then
\begin{equation}
\label{eq:thm-pointwise-normality-3}
\sqrt{n}\frac{\widehat g(x) - g(x)}{\|s(x)\|}\to_d N(0,1)
\end{equation}
\end{theorem}
\Cref{thm:pointwise-normality} shows that the estimator proposed in \Cref{sec:setup} 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.
\subsection{Uniform Convergence}
\label{subsec:uniform-inference}
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 \Cref{assm:second-stage-assumptions}.
\begin{assumption}[Uniform Limit Theory]
\label{assm:uniform-limit-theory}
Let \(\bar\eps_k = \sup_{1\leq j\leq k} |\eps_j|\), \(\alpha(x) := p^k(x)/\|p^k(x)\|\), and let
\[
\xi_k^L := \sup_{\substack{x,x'\in\calX \\ x\neq x'}} \frac{\|\alpha(x) - \alpha(x')\|}{\|x-x'\|}
.\]
Further for any integer \(s\) let \(\bar\sigma_k^s = \sup_{x\in\calX}\E[|\bar\eps_k|^s|X=x]\). For some \(m > 2\) assume
\begin{enumerate}[(i)]
\item The regression errors satisfy \(\sup_{x\in\calX} \E[\max_{1 \leq i \leq n} |\bar\eps_{k,i}|^m \mid X = x] \lesssim_P n^{1/m}\)
\item The basis functions are such that (a) \(\xi_k^{2m/(m-2)}\log k/n \lesssim 1\), (b) \((\bar\sigma_k^2\vee\bar\sigma_k^m)\log\xi_k^L \lesssim \log k\), and (c) \(\log\bar\sigma_k^m\xi_k \lesssim \log k\).
\end{enumerate}
\end{assumption}
As before, \Cref{assm:uniform-limit-theory} is very similar to its analogue in \citet{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)\).
\begin{theorem}[Strong Approximation by a Gaussian Process]
\label{thm:strong-approximation}
Assume that \Cref{cond:no-effect} holds and that Assumptions~\ref{assm:second-stage-assumptions}-\ref{assm:uniform-limit-theory} hold with \(m\geq 3\). In addition assume that (i) \(\bar R_{1n} = o_p(a_n^{-1})\) and (ii) \(a_n^6k^4\xi_k^2(\bar\sigma_k^3 + \ell_k^3c_k^2)^2\log^2 n/n\to 0\) where
\begin{align*}
\bar R_{1n} := \sqrt{\frac{\xi_k^2\log k}{n}}(n^{1/m}\sqrt{\log k} + \sqrt{k}\ell_k c_k) \andbox
\bar R_{2n} := \sqrt{\log k}\cdot\ell_kc_k
\end{align*}
Then so long as either the propensity score model or outcome regression model is correctly specified, for some \(\calN_k\sim N(0,I_k)\):
\begin{equation}
\label{eq:uniform-normality-beta}
\sqrt{n}\frac{\alpha(x)'(\widehat\beta-\beta)}{\|\alpha(x)'\Omega^{1/2}\|} =_d \frac{\alpha(x)'\Omega^{1/2}}{\|\alpha(x)'\Omega^{1/2}\|}N_k + o_p(a_n^{-1})\;\;\text{in }\ell^\infty(\calX)
\end{equation}
so that for \(s(x) := \Omega^{1/2}p^k(x)\)
\begin{equation}
\label{eq:uniform-normality-intermediate}
\sqrt{n}\frac{p^k(x)'(\widehat\beta-\beta)}{\|s(x)\|} =_d \frac{s(x)}{\|s(x)\|}N_k + o_p(a_n^{-1})\;\;\text{in }\ell^\infty(\calX)
\end{equation}
and if \(\sup_{x\in\calX} \sqrt{n}|r_k(x)|/\|s(x)\| = o(a_n^{-1})\), then
\begin{equation}
\label{eq:uniform-normality-ghat}
\sqrt{n}\frac{\widehat g(x) - g(x)}{\|s(x)\|} =_d \frac{s(x)'}{\|s(x)\|}\calN_k + o_p(a_n^{-1})\;\;\text{in }\ell^\infty(\calX)
\end{equation}
where in general we take \(\Omega = \tilde\Omega\) but if \(\bar R_{2n} = o_p(a_n^{-1})\) then we can set \(\Omega = \Omega_0\) where \(\tilde\Omega\) and \(\Omega_0\) are as in \eqref{eq:omega-definitions}.
\end{theorem}
\Cref{thm:strong-approximation} 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 \Cref{sec:setup}.
As noted by \citet{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.
\subsection{Matrix Estimation and Uniform Inference}
\label{subsec:uniform-bootstrap}
We establish that the estimator \(\widehat\Omega\) proposed in \eqref{eq:omega-hat-definitions} 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 \Cref{assm:second-stage-assumptions,assm:uniform-limit-theory} as well as the following condition limiting the first stage estimation error passed on to the variance estimator \(\widehat\Omega\).
\begin{condition}[Variance Estimation]
\label{cond:variance-estimation}
Let \(m > 2\) be as in \Cref{assm:uniform-limit-theory}. Then,
\begin{equation}
\label{eq:variance-condition}
\xi_{k,\infty}\max_{1\leq j\leq k}\E_n[p_j(X)^2(Y(\widehat\pi_j, \widehat m_j) - Y(\bar\pi_j, \bar m_j))^2] = o_p(k^{-2}n^{-1/m})
\end{equation}
\end{condition}
Via \Cref{thm:first-stage-convergence} we can establish \Cref{cond:variance-estimation} under \Cref{assm:logistic-model-convergence} as well as the additional sparsity bound\footnote{The sparsity bound \eqref{eq:variance-sparsity-bound} 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 \citet{SC-2020} and control alternate quantities with bounds that converge more quickly to zero.}
\begin{equation}
\label{eq:variance-sparsity-bound}
\frac{\xi_{k,\infty}^5s_k^2k^2\ln(d_z)}{n^{(m-1)/m}}
.\end{equation}
\begin{theorem}[Matrix Estimation]
\label{thm:matrix-estimation}
Suppose that Conditions~\ref{cond:no-effect} and \ref{cond:variance-estimation} and Assumptions~\ref{assm:second-stage-assumptions}-\ref{assm:uniform-limit-theory} hold. In addition, assume that \(\bar R_{1n} + \bar R_{2n} \lesssim (\log k)^{1/2}\). Then, so long as either the propensity score model or outcome regression model is correctly specified then for \(\widehat\Omega = \widehat{Q}^{-1}\widehat\Sigma \widehat{Q}^{-1}\):
\begin{align*}
\|\widehat\Omega - \Omega\| &\lesssim_P (v_n \vee \ell_kc_k)
\sqrt{\frac{\xi_k^2\log k}{n}} = o(1)
\end{align*}
\end{theorem}
\Cref{thm:matrix-estimation} establishes that pointwise inference based on the test statistic described in \Cref{sec:setup}, obtained by replacing \(\Omega\) in \Cref{thm:pointwise-normality} 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 \Cref{sec:setup}.
\begin{theorem}[Validity of Uniform Confidence Bands]
\label{thm:uniform-confidence-bands}
Suppose \Cref{cond:no-effect,cond:variance-estimation} are satisfied and \Crefrange{assm:second-stage-assumptions}{assm:uniform-limit-theory} hold with \(m\geq 4\). In addition suppose (i) \(R_{1n} + R_{2n} \lesssim \log^{1/2} n\), (ii) \(\xi_k\log^2 n /n^{1/2-1/m} = o(1)\), (iii) \(\sup_{x\in\calX} |r_k(x)|/\|p^k(x)\| = o(\log^{-1/2} n)\), and (iv) \(k^4\xi_k^2(1 + l_k^3r_k^3)^2\log^5 n/n=o(1)\). Then, so long as either the propensity score model or outcome regression model is satisfied
\[
\Pr\left(\sup_{x\in\calX} |\frac{\widehat g(x)-g(x)}{\widehat\sigma(x)}| \leq c^\star(1-\alpha)\right) = 1-\alpha + o(1)
.\]
As a result, uniform confidence intervals formed in \eqref{eq:confidence-bands} satisfy
\[
\Pr(g(x) \in [\underline i(x), \bar i(x)],\;\forall x\in\calX) = 1 - \alpha + o(1)
.\]
\end{theorem}
In conjunction with \Cref{thm:first-stage-convergence}, \Cref{thm:pointwise-normality} and \Cref{thm:matrix-estimation}, \Cref{thm:uniform-confidence-bands} shows the validity of the uniform inference procedure described in \Cref{sec:setup}.
\section{Estimation of the Conditional Average Treatment Effect}
\label{sec:cate-wrapup}
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 \eqref{eq:nuisance-parameter-functional-forms} 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 \eqref{eq:gamma-j-estimating-equation}-\eqref{eq:alpha-j-estimating-equation} however, and are instead given
\begin{align*}
\widehat\gamma_j^0 &:= \arg\min_\gamma\, \E_n[p_j(X)\{(1-D)e^{\gamma'Z} - D\gamma'Z\}] + \lambda_{\gamma,j}\|\gamma\|_1 \\
\widehat\alpha_j^0 &:= \arg\min_\alpha\, \E_n[p_j(Z)(1-D)e^{{\widehat{\gamma}_j^{0'}}Z}(Y - \alpha'Z)^2]/2 + \lambda_{\alpha,j}\|\alpha\|_1
\end{align*}
which under the natural analog of \Cref{assm:logistic-model-convergence} converge uniformly to population minimizers:
\begin{align*}
\bar\gamma_j^0 &:= \arg\min_\gamma\, \E[p_j(X)\{(1-D)e^{\gamma'Z} - D\gamma'Z\}] \\
\bar\alpha_j^0 &:= \arg\min_\alpha\, \E[p_j(Z)(1-D)e^{{\bar{\gamma}_j^{0'}}Z}(Y - \alpha'Z)^2]
\end{align*}
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 \Cref{sec:theory-overview}, 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 \Cref{sec:theory-overview,sec:first-stage}. 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 \eqref{eq:omega-definitions} but \(\Omega_1\) and \(\Omega_2\) are given
\begin{equation}
\label{eq:omega12}
\begin{split}
\Omega_1 &= Q^{-1}\E[\{p^k(x)\circ\eps_0^k\}\{p^k(x)\circ\eps_0^k\}']Q^{-1} \\
\Omega_2 &= Q^{-1}\E[\{p^k(x)\circ\eps^k\}\{p^k(x)\circ\eps_0^k\}']Q^{-1}
\end{split}
\end{equation}
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 \eqref{eq:omega-hat-definitions}.
\section{Simulation Study}
\label{sec:simulations}
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.
\subsection{Simulation Design}
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
\begin{align*}
\mathrm{Cov}(Z_{1,j},Z_{1,k}) = \E [Z_{1,j} Z_{1,k}] = 2^{-|j-k|}, \ \ \ 3\leq j,k \leq d_z.
\end{align*}
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 \texttt{sparsity} control the number of regressors in \(Z = (Z_1,X)\) entering the DGP.
\begin{enumerate}[label=(S\arabic*)]
\item \textit{Correct specification}: Generate \(D\) given \(Z\) from a Bernoulli distribution with \(\Pr (D=1 | Z ) = \{ 1+\exp(p_1 - X - 0.5X^2 - \gamma'Z_1) \}^{-1}\) and \(Y = D(1 + X + 0.5X^2 + \gamma'Z_1) + \epsilon.\)
\item \textit{Propensity score model correctly specified, but outcome regression model misspecified}: Generate \(D\) given \(Z\) as in (S1), but \( Y = D(1 + X + 0.5X^2 + \gamma'Z_1^\dagger) + \epsilon. \)
\item \textit{Propensity score model misspecified, but outcome regression model correctly specified}: Generate \(Y\) according to (S1), but generate \(D\) given \(Z\) from a Bernoulli distribution with \(\Pr (D=1 | Z ) = \{ 1+\exp(p_2 - X - 0.5X^2 + \gamma'Z_1^\dagger) \}^{-1}\).
\end{enumerate}
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), \texttt{sparsity}\(=6\); for (S2), \texttt{sparsity}\(=4\); and, for (S3), \texttt{sparsity}\(=5\). Results are reported for \(S=1,000\) repeated simulations.
\subsection{Estimators and Implementation}
To select the first stage penalty parameters, we implement the multiplier bootstrap procedure described in \Cref{subsec:additional}. The constants \(c_{\gamma,j}\) and \(c_{\alpha,j}\) in the pilot penalty parameters \eqref{eq:pilot-penalty} 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 \texttt{splines2} \citep{splines2-paper}, which uses the specification detailed in \cite{perperoglou2019review}. In the tables below, we refer to our method as \textit{MA-DML} (model assisted double machine learning).
We compare our proposed estimator and inference procedure to that of \citet{SC-2020}, which projects a single aIPW signal onto a growing series of basis terms. In implementing this \textit{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 \texttt{R} scripts implementing \textit{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), \)
\begin{align*}
\mathrm{IBias}^2 &= \int_0^1 (\Bar{g} (x) - g_0 (x))^2 dx, \\
\mathrm{IVar} &= S^{-1} \sum_{s=1}^S \int_0^1 (\widehat{g}_s (x) - \Bar{g} (x))^2 dx, \\
\mathrm{IMSE} &= S^{-1} \sum_{s=1}^S \int_0^1 (\widehat{g}_s (x) - g_0 (x))^2 dx. \\
\end{align*}
\subsection{Simulation Results}
Table \ref{tab:simulation} 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).
\begin{table}
\caption{Simulation study.} \label{tab:simulation}
\begin{center}
\begin{tabular}{ p{1cm} l c c c c c c c}
\hline \hline \\[-10pt]
DGP & Estimator & IBias$^2$ & IVar & IMSE & Cov90 & Cov95 & UCov90 & UCov95 \\[1pt]
& & (1) & (2) & (3) & (4) & (5) & (6) & (7) \\[1ex]
\hline
& & \multicolumn{7}{c}{K=3, n=500, $d_z$ = 100} \\ \cline{3-9}
\multirow{2}{1em}{(S1)} & DML & 0.04 & 0.31 & 0.35 & 0.92 & 0.96 & 1.00 & 1.00 \\
& MA-DML & $\sim$0.0 & 0.34 & 0.34 & 0.93 & 0.97 & 1.00 & 1.00 \\[7pt]
\multirow{2}{1em}{(S2)} & DML & 0.16 & 2.17 & 2.33 & 0.92 & 0.97 & 0.83 & 0.86 \\
& MA-DML & 0.03 & 2.12 & 2.15 & 0.90 & 0.94 & 0.88 & 0.91 \\ [7pt]
\multirow{2}{1em}{(S3)} & DML & 0.03 & 0.55 & 0.59 & 0.87 & 0.93 & 0.95 & 0.97 \\
& MA-DML & 0.01 & 0.79 & 0.80 & 0.91 & 0.95 & 0.99 & 0.99 \\[1ex] \cline{3-9}
& & \multicolumn{7}{c}{K=3, n=1000, $d_z$ = 100} \\ \cline{3-9}
\multirow{2}{1em}{(S1)} & DML & 0.12 & 0.20 & 0.32 & 0.83 & 0.90 & 0.96 & 0.96 \\
& MA-DML & 0.01 & 0.22 & 0.23 & 0.83 & 0.90 & 0.99 & 0.99 \\[7pt]
\multirow{2}{1em}{(S2)} & DML & 0.40 & 2.1 & 2.5 & 0.84 & 0.91 & 0.33 & 0.39 \\
& MA-DML & 0.19 & 2.07 & 2.26 & 0.83 & 0.89 & 0.50 & 0.55 \\ [7pt]
\multirow{2}{1em}{(S3)} & DML & 0.11 & 0.34 & 0.46 & 0.74 & 0.82 & 0.80 & 0.84 \\
& MA-DML & 0.01 & 0.53 & 0.54 & 0.84 & 0.89 & 0.89 & 0.91 \\
\hline\hline \\[-5pt]
\multicolumn{9}{c}{\parbox{12.2cm}{\footnotesize Note: DGP refers to the three various data generating processes introduced above. IBias$^2$, IVar, and IMSE refer to integrated squared bias, variance, and mean squared error, respectively. Cov90, Cov95, UCov90, and UCov95 refer to the coverage proportion of the 90\% and 95\% pointwise and uniform confidence intervals across simulations. $K$ refers to the number of series terms, $N$ to the sample size, and $d_z$ to the dimensionality of the random variable $Z_1.$}} \\
\end{tabular}
\end{center}
\end{table}
For pointwise and uniform coverage under correct specification regime (S1), \textit{MA-DML} has some slight improvements. Under misspecification DGPs (S2) and (S3), the pointwise coverage of \textit{MA-DML} is closer to the targets except in the $N=1000$ and (S2) case where it slightly underperforms. However, \textit{MA-DML} has a notable improvement over \textit{DML} in the (S3) case when $N=1000.$ Similarly, \textit{MA-DML} outperforms \textit{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, \textit{MA-DML} has a lower MSE than \textit{DML} where regardless of sample size \textit{MA-DML} underperforms in (S3). Notably, it does appear \textit{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 \textit{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 \textit{MA-DML} is lower across the various DGPs compared to \textit{DML}. Following the discussion in \Cref{sec:theory-overview} 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 \citet{SC-2020} benchmark method, whose work we rely on and were inspired by.
\section{Empirical Application}
\label{sec:empirical}
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 \citet{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 \Cref{sec:empirical-details}.
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 \citet{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 \citep{Zimmert_Lechner_CATE,Abreya_2006,Lee_Ryo_Wang_2016}.
\subsection{Empirical Results}
\Cref{fig:emp1} 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 \Cref{sec:first-stage}. 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 \eqref{eq:final-penalty-parameters} 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 \texttt{R} \citep{splines2-paper}.
\begin{figure}[htp!]
\centering
\includegraphics[width=\linewidth]{figures/empirical/K4_D2_alpha0.01,0.05,0.1_5.png}
\caption{CATE of maternal smoking estimated using model assisted estimating equations (left) and standard MLE/OLS estimating equations (right). Top row uses the 99\textsuperscript{th} quantile of the bootstrap distribution to select the penalty parameters, second row uses 95\textsuperscript{th} quantile, and final row uses the 90\textsuperscript{th} quantile. Second stage is computed using b-splines of the second degree with 3 knots. 95\% pointwise confidence intervals are displayed in blue short dashes and 95\% uniform confidence bands are displayed in long red dashes.}
\label{fig:emp1}
\end{figure}
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 \eqref{eq:gamma-j-estimating-equation}-\eqref{eq:alpha-j-estimating-equation} 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 \citet{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 \citet{Zimmert_Lechner_CATE}, is the existence of the uniform confidence bands displayed. Reassuringly, the estimates of \citet{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 \Cref{fig:empD1}. 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.
\begin{figure}[htpb]
\centering
\includegraphics[width=\linewidth]{figures/empirical/K4_D1.png}
\caption{CATE of maternal smoking estimated using model assisted estimating equations (left) and standard MLE/OLS estimating equations (right). Top row uses the 99\textsuperscript{th} quantile of the bootstrap distribution to select the penalty parameters, second row uses 95\textsuperscript{th} quantile, and final row uses the 90\textsuperscript{th} quantile. Second stage is computed using b-splines of the first degree with 3 knots. 95\% pointwise confidence intervals are displayed in blue short dashes and 95\% uniform confidence bands are displayed in long red dashes.}
\label{fig:empD1}
\end{figure}
Finally, \Cref{tab:implied-ate} reports the smoothed average treatment effect estimates taken from averaging the model assisted CATE estimates from \Cref{fig:emp1} across observations. Again, these estimates are generally in line with prior work
\begin{table}[htpb]
\centering
\caption{Smoothed Model Assisted ATE Estimates}
\vspace{0.1em}
\label{tab:implied-ate}
\begin{tabular}{c|ccc}
\hline
\hline
Bootstrap Penalty Qt.& 99\textsuperscript{th} & 95\textsuperscript{th} & 90\textsuperscript{th} \\
\hline
Implied ATE & -295.221 & -292.9086 & -453.2242
\end{tabular}
\end{table}
\section{Conclusion}
\label{sec:conclusion}
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.
\newpage
\singlespacing
\bibliography{bibtex/mlci.bib}