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.
61,935 characters
Doubly Robust Uniform Confidence Band for the Conditional Average Treatment Effect Function
\doublespacing
\begin{abstract}
In this paper, we propose a doubly robust method to estimate the heterogeneity of the average treatment effect with respect to observed covariates of interest.
We consider a situation where a large number of covariates are needed for identifying the average treatment effect but the covariates of interest for analyzing heterogeneity are of much lower dimension.
Our proposed estimator is doubly robust and avoids the curse of dimensionality.
We propose a uniform confidence band that is easy to compute, and we illustrate its usefulness via Monte Carlo experiments and an application to the effects of smoking on birth weights.
\bigskip
\noindent
Keywords: average treatment effect conditional on covariates, uniform confidence band, double robustness, Gaussian approximation.
\bigskip
\noindent
JEL classification codes: C14, C21
\end{abstract}
\maketitle
\clearpage
\section{Introduction}
In this paper, we propose a doubly robust method to estimate the heterogeneity of the average treatment effect
with respect to observed covariates of interest. To describe our methodology,
we consider the potential outcome framework.
Let $Y_{1}$ and $Y_{0}$ be potential individual outcomes in two states, with
treatment and without treatment, respectively. For each individual, the observed outcome $
Y $ is $Y=DY_{1}+(1-D)Y_{0}$, where $D$ denotes an indicator
variable for the treatment, with $D=0$ if an individual is not treated and $
D=1$ if an individual is treated. We assume that independent and identically
distributed observations $\{(Y_{i},D_{i},\mathbf{Z}_{i}):i=1,\ldots ,n\}$ of $(Y,D,\mathbf{Z})$
are available, where $\mathbf{Z} \in \mathbb{R}^{p}$ denotes a $p$-dimensional vector of covariates.
Suppose that a researcher is interested in evaluating the average treatment effect conditional on only a
subset of covariates $\mathbf{X}$,
which is of a substantially lower dimension than $\mathbf{Z}$,
where $\mathbf{Z}\equiv (\mathbf{X}^{\top },\mathbf{V}^{\top })^{\top }\in \mathbb{R}^{d}\times
\mathbb{R}^{m}$, $p \equiv d+m$.
That is, we are interested in a case where $d \ll p$.
The main object of interest in this paper is the conditional average treatment effect function (CATEF); namely:
\begin{align}\label{def:cater}
g(\mathbf{x}) \equiv\mathbb{E}[Y_{1}-Y_{0}|\mathbf{X}=\mathbf{x}].
\end{align}
When $d \geq 3$, it is difficult to plot $g(\mathbf{x})$, not to mention low precision due to the curse of dimensionality.
Hence, for practical reasons, we focus on the case that $d=1$ or $d=2$, while $p$ is often of a much higher dimension.
To achieve identification of the CATEF, we
assume that $Y_1$ and $Y_0$ are independent of
$D$ conditional on $\mathbf{Z}$ (known as the unconfoundedness assumption):
\begin{align}\label{unconfound}
(Y_1, Y_0) \perp D |\mathbf{Z},
\end{align}
where $\perp$ denotes the independence.
For \eqref{unconfound} to be plausible in applications, applied researchers
tend to consider a large number of covariates $\mathbf{Z}$.
Note that in our setup, the treatment may be confounded in the sense that the treatment assignment may not be independent of the potential outcome variables given $\mathbf{X}$ only.
To satisfy the unconfoundedness condition, a much larger set of conditioning variables $\mathbf{Z}$ needs to be employed.
Different roles of covariates between $\mathbf{X}$ and $\mathbf{V}$
are noted in the recent literature.
For example, \cite{Ogburn-et-al-JRSSB} consider a similar issue in the context of
the local average treatment effect (LATE) of \cite{Imbens:Angrist:94}.
\cite{Ogburn-et-al-JRSSB} emphasize that
conditioning on a large number of covariates $\mathbf{Z}$ may be required to make it plausible that
the binary instrument is valid.
In their empirical example, \cite{Ogburn-et-al-JRSSB} revisit the analyses of
\cite{Poterba1995} and \cite{Abadie2003} to examine whether participation in 401(k) pension plans increases household savings.
In their example, the vector of covariates $\mathbf{Z}$ for the identifying assumption consists of
income, age, marital status, and family size, whereas the variable of interest $\mathbf{X}$ is income.
\cite{Abrevaya-et-al-JBES} also consider the case of investigating the effect of smoking during pregnancy on birth weights.
They are interested in estimating \eqref{def:cater} with $\mathbf{X}$ being the age of mother; however, as noted in \cite{Abrevaya-et-al-JBES}, it is unlikely that conditioning only on the age of mother would achieve the unconfoundedness assumption with nonexperimental data. As a result,
it is necessary to consider a high-dimensional $\mathbf{Z}$, including the age of the mother.
The fact that a high-dimensional $\mathbf{Z}$ needs to be employed for \eqref{unconfound} to be plausible in an application makes a fully
nonparametric estimation approach impractical because of the curse of
dimensionality. For example, the propensity score is not nonparametrically estimable
in moderately sized samples, if the dimension of $\mathbf{Z}$ is high.
One obvious alternative is to use a parametric model for the propensity score; however, it may lead to misleading results if the parametric model is misspecified.
With the aim of providing a practical method and, at the same time, reducing sensitivity to model misspecification,
we propose to use a doubly robust method based on parametric regression and propensity score models.
Our estimator of the CATEF is doubly robust in the sense that it is consistent when at least one of the regression model and the propensity score model is correctly specified.
Specifically, we first estimate CATEF$(\mathbf{Z})$ using a doubly robust procedure:
we estimate a parametric regression model of the outcome on $\mathbf{Z}$ for each treatment status and a parametric model for the probability of selecting into the treatment given $\mathbf{Z}$; we then
combine the parametric estimation results in a doubly robust fashion to construct an estimate of CATEF$(\mathbf{Z})$.
We then obtain an estimate of CATEF$(\mathbf{X})$ by adopting the local linear smoothing of CATEF$(\mathbf{Z})$.
As a result, we avoid high-dimensional smoothing with respect to
$\mathbf{Z}$ but mitigate the problem of misspecification by both the doubly robust estimation and low-dimensional nonparametric smoothing with respect to $\mathbf{X}$.
We emphasize that we are willing to assume parametric specifications for the propensity score and
regression models as functions of $\mathbf{Z}$ to avoid the curse of dimensionality, but not for CATEF$(\mathbf{X})$.
One may consider parametric estimation of CATEF$(\mathbf{X})$, as \cite{Ogburn-et-al-JRSSB} estimate their LATE parameter
using least squares approximations. However, note that
even if the parametric specification of CATEF$(\mathbf{Z})$ is correct, the resulting specification of CATEF$(\mathbf{X})$
may not be correctly specified since, for example, $\mathbb{E}[ \mathbf{Z}| \mathbf{X}]$ is possibly highly nonlinear.
To avoid this misspecification, we estimate CATEF$(\mathbf{X})$ nonparametrically.
Because the CATEF is a functional parameter, as a tool of inference, we propose to use a uniform confidence band for the CATEF.
Our construction of the uniform confidence band is based on some analytic approximation
of the supremum of a Gaussian process using arguments built on \cite{Piterbarg:96}, combined with a Gaussian approximation result of
\cite{CCK:2014} and an empirical process result of \cite{Ghosal/Sen/vanderVaart:00}.
Our method is simple to implement
and does not rely on resampling techniques.
This paper contributes to the literature on doubly robust estimation by demonstrating that the doubly robust procedures are useful for estimating the CATEF.
In this paper, we focus on the so-called augmented inverse probability weighting estimator that was originally proposed by \cite{RobinsRotnitzkyZhao94} for the estimation of the mean \citep[see also][]{RobinsRotnitzky95, ScharfsteinRotnitzkyRobins99}.
Their estimator appears to be the first estimator to be recognized as being doubly robust.
Since then, many other alternative doubly robust estimators have been proposed in the literature.
For example, the inverse probability weighting regression adjustment estimator \citep{KangSchafer07, Wooldridge07, Wooldridge10}
is widely known and has been implemented in statistical software packages.
See the introduction of \cite{Tan10} for a comprehensive summary of other doubly robust estimators.
Doubly robust estimators have been advocated for use in many different areas of application: See, for example, \cite{LuncefordDavidian04} for medicine, \cite{GlynnQuinn10} for political science, \cite{Wooldridge10} for economics, and \cite{SchaferKang08} for psychology.
There are also doubly robust estimators available for different settings including instrumental variables estimation \citep{Tan06, OkuiSmallTanRobins12}
and estimation under multivalued treatments \citep{DeryaUysal:15}.
It would not be difficult to extend our method to allow these other doubly robust estimators and to consider different settings.
However, to keep the analysis simple, in this paper, we focus on the augmented inverse probability weighting
estimator of the CATEF.
The CATEF is mathematically equivalent to ``$V$-adjusted variable importance'' of \cite{vanderLaan06}, who proposes it as a measure of variable importance in prediction.
\cite{vanderLaan06} proposes a doubly robust estimator of $V$-adjusted variable importance.
Contrary to ours, he considers the projection of the $V$-adjusted variable importance
on a parametric working model and
does not consider a nonparametric estimation.
Moreover, a uniform confidence band is not examined in \cite{vanderLaan06}.
In a recent paper, \cite{Abrevaya-et-al-JBES} consider the estimation of the CATEF\footnote{Our paper is independent of \cite{Abrevaya-et-al-JBES} and it is started without knowing their work.}; however, there are two main differences of this paper relative to \cite{Abrevaya-et-al-JBES}.
First, we propose the doubly robust procedure to estimate the CATEF.
\cite{Abrevaya-et-al-JBES} consider the inverse probability weighting estimator.
The inverse probability weighting estimator suffers from model misspecification when the propensity score model is misspecified
and from the curse of dimensionality when it is estimated nonparametrically.
Second, we present a method to construct a uniform confidence band, whereas
\cite{Abrevaya-et-al-JBES} only provide a pointwise confidence interval.
The remainder of the paper is organized as follows.
Section \ref{sec:method} presents the doubly robust estimation method,
Section \ref{sec:simple:case} gives an informal description of how to construct a two-sided, symmetric uniform confidence band when the dimension of $\mathbf{X}$ is one,
and Section \ref{sec:general:case} deals with a general case and provides formal theoretical results. In Section \ref{sec:MC}, the results of Monte Carlo simulations demonstrate that in finite samples, our doubly robust estimator works well, and the proposed confidence band has desirable coverage properties.
Section \ref{sec:EA} gives an empirical application, and Section \ref{sec:C} concludes.
The proofs are contained in Appendix \ref{sec:A}.
\section{Doubly Robust Estimation of the Average Treatment Effect Conditional on
Covariates of Interest}\label{sec:method}
In this section, a doubly robust method for estimating the CATEF is proposed.
We first estimate the CATEF for all the covariates using a doubly robust method.
We then obtain the CATEF for the covariates of interest using a nonparametric approach.
Define:
\begin{align*}
\pi (\mathbf{z}) & \equiv \mathbb{E} \left[ D|\mathbf{Z}=\mathbf{z}\right], \\
\mu _{j}(\mathbf{z}) & \equiv \mathbb{E} \left[ Y|\mathbf{Z}=\mathbf{z},\ D=j\right] \ \text{for }j=0,1,
\end{align*}
where $\pi (\mathbf{z})$ is the propensity score and $\mu_j (\mathbf{z})$ for $j=0,1$ are called regression functions. Note that $\mu_j (\mathbf{z}) = E( Y_j | \mathbf{Z}= \mathbf{z})$ for $j=0,1$ under unconfoundedness.
Let $\pi (\mathbf{z},\beta )\ $and $\mu _{j}(\mathbf{z},\alpha _{j})$ for $j=0,1$ denote
parametric models of $\pi (\mathbf{z})$ and $\mu _{j}(\mathbf{z})$,
respectively.\footnote{$\mu _{j}(\mathbf{z},\alpha _{j})$ may also be called ``marginal structural models'' of \cite{Robins00}.}
A doubly robust procedure requires that either $\pi (\mathbf{z})$ or $\mu _{j}(\mathbf{z})$ for $j=0,1$
should be correctly specified, thereby allowing for misspecification in $\pi (\mathbf{z})$ or in $\mu _{j}(\mathbf{z})$.
Let $\theta_0 \equiv ( \alpha_{10}^\top, \alpha_{00}^\top,\beta_0^\top)^\top$ denote the vector of true or pseudo-true parameter values that optimize some criterion functions.
We consider the augmented inverse probability weighting approach.
Let:
\begin{eqnarray*}
\psi _{1}(\mathbf{W},\alpha _{1},\beta ) & \equiv&\frac{DY}{\pi (\mathbf{Z},\beta )}-\frac{D-\pi
(\mathbf{Z},\beta )}{\pi (\mathbf{Z},\beta )}\mu _{1}(\mathbf{Z},\alpha _{1}), \\
\psi _{0}(\mathbf{W},\alpha _{0},\beta ) & \equiv&\frac{\left( 1-D\right) Y}{1-\pi (\mathbf{Z},\beta
)}+\frac{D-\pi (\mathbf{Z},\beta )}{1-\pi (\mathbf{Z},\beta )}\mu _{0}(\mathbf{Z},\alpha _{0}), \\
\psi (\mathbf{W},\theta ) & \equiv&\psi _{1}(\mathbf{W},\alpha _{1},\beta )-\psi _{0}(\mathbf{W},\alpha
_{0},\beta ),
\end{eqnarray*}
where $\mathbf{W} \equiv (Y,\mathbf{Z}^{\top })^{\top }\ $and $\theta \equiv (\alpha _{1}^{\top
},\alpha _{0}^{\top },\beta ^{\top })^{\top }.$
The first terms in $\psi_{1}(\mathbf{W},\alpha _{1},\beta ) $ and $\psi_{0}(\mathbf{W},\alpha _{0},\beta ) $
correspond to inverse probability weighting.
The second terms are augmented terms that make the procedure doubly robust.
The following lemma gives
regularity conditions under which
$g(\mathbf{x})$ is identified.
\begin{lem}[Identification of the CATEF]\label{iden-lem}
Assume that \eqref{unconfound} holds and $0 < \pi(\mathbf{Z},\beta_0 ) < 1$ almost surely.
Suppose that either $\beta_0$ satisfies $\mathbb{E}\left[D|\mathbf{Z}\right] = \pi
(\mathbf{Z},\beta_0 )$ almost surely or $\alpha_{10}$ and $\alpha_{00}$ satisfy
$\mathbb{E}\left[Y_1|\mathbf{Z}\right] = \mu _{1}(\mathbf{Z},\alpha _{10})$
and $\mathbb{E}\left[Y_0|\mathbf{Z}\right] = \mu _{0}(\mathbf{Z},\alpha _{00})$ almost surely.
Then:
\begin{align*}
g(\mathbf{x}) = \mathbb{E}\left[ \psi (\mathbf{W},\theta_0 )|\mathbf{X}=\mathbf{x} \right].\end{align*}
\end{lem}
Lemma \ref{iden-lem} suggests that one may estimate
$g(\mathbf{x})$ by running the nonparametric regression of $\psi (\mathbf{W}, \hat \theta)$
on $\mathbf{X}_i$, where $\hat \theta$ is a consistent parametric estimator of $\theta_0$.
Moreover, this lemma implies that the CATEF can be identified through $\psi (\mathbf{W}, \theta_0)$ if either the regression models ($\mu_1 (\mathbf{z}, \alpha_1)$ and $\mu_0 (\mathbf{z}, \alpha_0)$) or the propensity score model ($\pi (\mathbf{z},\beta)$) is correctly specified (or both).
That is, even if $\mu_1 (\mathbf{z}, \alpha_1)$ and $\mu_0 (\mathbf{z}, \alpha_0)$ do not represent the true conditional expectation functions, provided that $\pi (\mathbf{z}, \beta)$ is correct, the CATEF is identified. Similarly, even if $\pi (\mathbf{z}, \beta)$ is misspecified,
provided that $\mu_1 (\mathbf{z}, \alpha_1)$ and $\mu_0 (\mathbf{z}, \alpha_0)$ are correct,
the CATEF is identified.
\begin{rem}
In this paper, we focus on cases in which $\mathbf{X}$ is continuous. When $\mathbf{X}$ is discrete, the CATEF can be estimated by the sample average of $\psi (\mathbf{W}, \hat \theta)$ using the sub-sample for each possible value of $\mathbf{X}$ and an estimator $\hat \theta$ of $\theta_0$. Moreover, constructing a confidence band is standard when $\mathbf{X}$ takes a finite number of values.
\end{rem}
\subsection{Parametric Estimation of $\theta$}
For concreteness, we consider the following estimation procedure for $\theta_0$.
However, how $\theta_0$ is estimated does not alter our results provided that the rate of convergence is sufficiently fast so that Assumption \ref{main-assumption}(7) given below is satisfied.
For each $j = 0,1$, we estimate $\alpha_j$ by least squares:
\begin{align}\label{ols-est}
\hat \alpha_j \equiv \operatorname*{argmin}_{\alpha_j}
\sum_{i=1}^n D_i^j (1-D_i)^{1-j} [Y_i - \mu_j (\mathbf{Z}_i, \alpha_j)]^2.
\end{align}
We estimate $\beta$ by maximum likelihood (e.g., probit or logit):
\begin{align}
\hat \beta \equiv \operatorname*{argmax}_{\beta} \sum_{i=1}^n \left(
D_i \log \pi(\mathbf{Z}_i, \beta) + (1-D_i) \log (1- \pi (\mathbf{Z}_i,\beta))\right).
\end{align}
\begin{rem}
When the dimension of $\mathbf{Z}$ is not too high, an alternative to parametric estimation of
$\psi (\mathbf{W},\theta_0 )$ is to estimate its nonparametric counterpart via local polynomial estimators
as in \cite{Rothe/Firpo:16}. However, this would not work when the dimension of $\mathbf{Z}$ is sufficiently high (see related remarks in \cite{Rothe/Firpo:16}).
The latter is the case we focus on in the paper.
\end{rem}
\subsection{Local Linear Estimation of $g$}
We consider a
local linear estimator of $g(\mathbf{x})$.
Assume that $g(\mathbf{x})$ is twice
continuously differentiable.
For each $\mathbf{x}=(x_{1},\ldots ,x_{d})$, the local linear estimator of $g(\mathbf{x})$ can be obtained by
minimizing:
\begin{equation*}
S_{n}(\gamma )\equiv \sum_{i=1}^{n}\left[ \psi (\mathbf{W}_{i},\hat \theta
)- \gamma_0 - \gamma_1 ^{\top }\left(\mathbf{X}_i - \mathbf{x} \right) \right]^2
\mathbf{K} \left( \frac{\mathbf{X}_i - \mathbf{x}}{h_n} \right)
\end{equation*}
with respect to $\gamma \equiv (\gamma_0, \gamma_1^\top) ^\top \in \mathbb{R}^{d+1}$, where $\mathbf{K}(\cdot)$ is a kernel function on $\mathbb{R}^d$
and $h_n$ is a sequence of bandwidths. More specifically, let $
\hat{g}( \mathbf{x})=\mathbf{e}_{1}^{\top }\hat{\gamma}(\mathbf{x})$, where $
\hat{\gamma}(\mathbf{x})\equiv \arg \min_{\gamma \in \mathbb{R}
^{d+1}}S_{n}(\gamma )$ and $\mathbf{e}_{1}$ is a column vector
whose first entry is one, and the rest are zero.
\subsection{Effect of First Stage Estimation}
In our setting, we can carry out inference as if $\theta_0$ were known.
This result would not be a surprise given that our first-stage estimation is parametric and our second-stage estimation is nonparametric: the rate of the convergence in the first-stage estimation is faster than that of the second stage.
This feature of no first-order effect of the first-stage estimation in the second stage turns out to be more general than our setup.
It is indeed closely related to doubly robustness.
If we model $g(\mathbf{x})$ parametrically or more generally approximate $g(\mathbf{x})$ by linear projection, it can be estimated by running an OLS of $\psi (\mathbf{W},\hat \theta )$ on $\mathbf{X}$.
Because of the built-in feature of double robustness, it can be shown that the limiting distribution of
the OLS estimator of $\psi (\mathbf{W},\hat \theta )$ on $\mathbf{X}$ is equivalent to that of the infeasible
OLS estimator of $\psi (\mathbf{W}, \theta_0 )$ on $\mathbf{X}$.
Furthermore, even if we estimate $\pi(\cdot)$ and $\mu_j(\cdot)$ $(j=0,1)$ nonparametrically when the dimension of $\mathbf{Z}$ is moderate, there will be no estimation effect from the first stage as well.
For example, see \cite{chen2008}, \cite{Rothe/Firpo:16} and \cite{CEIN2016} among others for related results.
\section{An Informal Description of a Uniform Confidence Band}\label{sec:simple:case}
In this section, we provide an informal description of how to construct a two-sided, symmetric uniform confidence band.
For simplicity, we focus on the leading case where $d=1$.
Let $\mathcal{I} \equiv [a,b]$ denote an interval of interest for which we build a uniform confidence band. Assume that $\mathcal{I}$ is a subset of the support of $X$.
We use nonbold $x$ to mean that $x$ is one-dimensional.
\begin{algorithm}
Carry out the following steps to construct a $(1-\alpha)$ uniform confidence band.
\begin{enumerate}
\item Obtain $\hat{g}(x)$ using a local linear estimator with
a bandwidth $h_n$ such that:
\begin{align*}
h_n = \widehat{h} \times n^{1/5} \times n^{-2/7},
\end{align*}
where $\widehat{h}$ is a commonly used optimal bandwidth in the literature
(for example, the plug-in method of \cite{Ruppert:Sheather:Wand:1995} which is explained in Appendix \ref{ap-bandwidth}).
We use the Gaussian kernel in our simulations and empirical application.
\item Obtain the pointwise standard error $\hat{s}(x)/(n h_n)^{1/2}$ of $\hat{g}(x)$
by constructing a feasible version of the asymptotic standard error formula:
\begin{align}\label{se-form}
\frac{\hat s(x)}{(nh_n)^{1/2}} \equiv \left\{ [n h_n \hat f_X(x)]^{-1} \int K^2(u) du \, \hat \sigma^2(x) \right\}^{1/2},
\end{align}
where $\hat{f}_X$ is the kernel
density estimator:
\begin{align*}
\hat{f}_X (x) = \frac{1}{nh_n} \sum_{i=1}^n K \left(\frac{X_i - x}{h_n} \right),
\end{align*}
and $\hat{\sigma}^2(x)$ is the conditional variance function estimator:
\begin{align}\label{sigma-x-est-dim1}
\hat{\sigma}^2 (x)
= \frac{1}{(n-\dim (\theta) ) h_n} \sum_{i=1}^N \frac{\hat U_i^2}{\hat{f}_X (x)} K \left( \frac{X_i - x}{h_n}\right).
\end{align}
Here, $\hat U_i = \psi ( \mathbf{W}_i, \hat \theta) - \hat g(X_i)$
and $\dim(\theta)$ is the dimension of $\theta$.
\item To compute a critical value $c(1-\alpha)$, define:
\begin{equation*}
\lambda \equiv \frac{ - \int K(u) K^{\prime \prime }(u) du }{\int K^2(u) du }.
\end{equation*}
Note that $\lambda =0.5$ if $K (\cdot)$ is the Gaussian kernel.\footnote{Note that $\lambda= 1.98$ for the biweight kernel and $\lambda= 2.5$ for the Epanechnikov kernel.}
Let:
\begin{equation*}
a_n \equiv a_{n}(\mathcal{I})= \left( 2\log (h_{n}^{-1}(b-a))+2\log \frac{\lambda ^{1/2}}{2\pi }
\right) ^{1/2}.
\end{equation*}
Now set the critical value for the two-sided symmetric uniform confidence band by:
\begin{equation*}
c(1-\alpha) \equiv \left( a_n^2 - 2 \log \{ \log [(1-\alpha )^{-1/2}] \} \right)^{1/2}.
\end{equation*}
\item For each $x \in \mathcal{I}$, we set the two-sided symmetric confidence band:
\begin{align*}
\hat{g}(x) - c(1-\alpha) \frac{\hat{s}(x)}{\sqrt{nh_n}}
\leq
g(x)
\leq
\hat{g}(x) + c(1-\alpha) \frac{\hat{s}(x)}{\sqrt{nh_n}}.
\end{align*}
\end{enumerate}
\end{algorithm}
\bigskip
We make some remarks on the proposed algorithm.
In step (1),
the factor $n^{1/5} \times n^{-2/7}$ is multiplied in the definition of $h_n$ to ensure that the bias is asymptotically negligible by undersmoothing.
In step (2), one can estimate $f_X$ and $\sigma^2(x) \equiv \text{Var}\left[ \psi (\mathbf{W},\theta_0 )|X=x \right]$ using the standard kernel density and regression estimators with the same kernel function $K(\cdot)$ and the same bandwidth $h$ and also with an estimator of $\theta_0$.
In step (3), we may restrict the bandwidth such that $h_{n} \ll (b-a)$ (which is satisfied asymptotically), thereby imposing the condition that $\log (h_{n}^{-1}(b-a))$ is positive.
The critical value proposed in step (3) is strictly positive if $\alpha$ is not too close to one or if $n$ is large enough.
\begin{rem}
It is straightforward to modify the algorithm above to construct one-sided symmetric confidence bands.
Define a new critical value by
\begin{equation*}
c_{\text{one-sided}}(1-\alpha) \equiv \left( a_n^2 - 2 \log \{ \log [(1-\alpha )^{-1}] \} \right)^{1/2}.
\end{equation*}
Then, for each $x \in \mathcal{I}$, we set the one-sided symmetric confidence bands:
\begin{align*}
\hat{g}(x) - c_{\text{one-sided}}(1-\alpha) \frac{\hat{s}(x)}{\sqrt{nh_n}}
\leq g(x),
\end{align*}
or
\begin{align*}
g(x)
\leq
\hat{g}(x) + c_{\text{one-sided}}(1-\alpha) \frac{\hat{s}(x)}{\sqrt{nh_n}}.
\end{align*}
\end{rem}
\begin{rem}{\label{rem-algorithm-multivariate}}
When $x$ is more than one dimension, the algorithm may be revised as follows.
Obviously, we need to use multivariate kernels and pointwise standard errors should be adjusted because the rate of convergence becomes $nh_n^d$.
The value of $\lambda$ stays the same when we use a product kernel. For example, if $\mathbf{K}$ is the product Gaussian kernel, then $\lambda =0.5$.
The formulas of $a_n$ and $c(1-\alpha)$ need to be changed.
$a_n$ is the largest solution to the following equation:
\begin{align*}
\text{mes}(\mathcal{I}) {h_n}^{-d} \lambda^{d/2} (2\pi)^{-(d+1)/2 }
a_n^{d-1} \exp(-a_n^2/2) = 1,
\end{align*}
where $\text{mes}(\mathcal{I})$ is the Lebesgue measure of $\mathcal{I}$.
When $d =2$, the critical value has the form $c(1-\alpha) \equiv a_n + c/a_n$, where $c$ is the smallest value that satisfies
\begin{align*}
\exp \left( -2 e^{-c -c^2 / 2 a_n^2}\right)\left( 1+ \frac{c}{a_n^2}\right) \ge 1-\alpha.
\end{align*}
When $d=3$, we have that $c(1-\alpha) \equiv a_n + c/a_n$, where $c$ is the smallest value that satisfies
\begin{align*}
\exp \left( -2 e^{-c -c^2 / 2 a_n^2}\right)\left( \left( 1+ \frac{c}{a_n^2}\right)^2 - 2 \frac{1}{a_n^2} \right) \ge 1- \alpha.
\end{align*}
We note that in this paper, we assume that $d <4$ (see Assumption \ref{main-assumption}).
\end{rem}
\begin{rem}\label{rem-gumbel}
We may compare our proposal with the critical value based on
the $(1-\alpha)$ quantile of the Gumbel distribution, which is given by:
\begin{equation*}
c_{\infty}(1-\alpha) \equiv a_{n} + \frac{-\log \{ \log [(1-\alpha)^{-1/2}] \}}{a_{n}}.
\end{equation*}
Note that:
$$
c_{\infty}(1-\alpha) - c(1-\alpha) = \left[\frac{-\log \{ \log [(1-\alpha)^{-1/2}] \}}{a_{n}}\right]^2,
$$
which is strictly positive for small $\alpha$ but converges to zero as $a_n$ diverges.
Hence, we expect that in finite samples, the confidence band based on $c_{\infty}(1-\alpha)$ is too wide and
has a higher coverage probability than the nominal level.
It is shown in the next section that the critical value based on the Gumbel distribution is accurate only up to the logarithmic rate, where our proposed critical value is precise in a polynomial rate.
This is because our proposal uses a higher-order expansion of \cite{Piterbarg:96},
whose approximation error is of a polynomial rate.
See Theorem \ref{main-uniform} in Section \ref{sec:general:case} for details.
\end{rem}
\begin{rem}
Our construction of critical values is based on a simple analytic method that is easy to compute.
Alternatively, one may rely on bootstrap methods to compute critical values for the uniform confidence band.
For example, see \cite{claeskens2003}
for smoothed bootstrap confidence
bands and \cite{CLR}
for multiplier bootstrap confidence bands.
\cite{chernozhukov2013}
show that in general settings including high dimensional models,
Gaussian multiplier bootstrap methods yield critical values for which the approximation error decreases polynomially in the sample size. Roughly speaking, both our simple analytic correction and multiplier bootstrap methods
yield critical values that are accurate at polynomial rates.
A refined theoretical analysis is necessary to determine which type of the critical value is better asymptotically.
\end{rem}
\begin{rem}
The proposed confidence band can be used to test whether the CATEF is constant.
Suppose that our null hypothesis is that $g(\mathbf{x})$ is constant in $ \mathcal{I}$.
This null hypothesis can be written as $g(\mathbf{x}) = g_{\mathcal{I}}$, where $g_{\mathcal{I}}
= \mathbb{E}[ g(\mathbf{x}) | \mathbf{x} \in \mathcal{I}]$.
Since $g_{\mathcal{I}}$ can be estimated at the parametric ($\sqrt{n}\,$) rate and the estimator thus converges faster than $\hat g(\mathbf{x})$, we can ignore the estimation error for $g_{\mathcal{I}}$.
We reject the constancy of $g(\mathbf{x})$, if the confidence band does not include the estimate of $g_{\mathcal{I}}$ for some $\mathbf{x} \in \mathcal{I}$.
\end{rem}
\section{Asymptotic Theory}\label{sec:general:case}
In this section, we establish asymptotic theory.
Let $U \equiv \psi (\mathbf{W},\theta_0 ) - g(\mathbf{X})$ and let
$U_i \equiv \psi (\mathbf{W}_i,\theta_0 ) - g(\mathbf{X}_i)$ for $i=1,\ldots,n$.
Let $\hat s^2 (\mathbf{x})$ be the estimator of the asymptotic variance of $\hat g(\mathbf{x})$.
Let $s_n^2(\mathbf{x})$ denote the population version of the asymptotic variance of the estimator:
\begin{align*}
s_n^2(\mathbf{x})
\equiv \frac{1}{h_n^{d}} \mathbb{E}\left[ \frac{U^2}{f^2_\mathbf{X}(\mathbf{x})} \mathbf{K}^2 \left( \frac{\mathbf{X} - \mathbf{x}}{h_n} \right) \right].
\end{align*}
Assume that the $d$-dimensional kernel function is the product of $d$ univariate kernel functions.\
That is, $\mathbf{K}(\mathbf{s}) = \prod_{j=1}^d K(s_j)$, where
$\mathbf{s} \equiv (s_1,\ldots,s_d)$ is a $d$-dimensional vector
and $K(\cdot)$ is a kernel function on $\mathbb{R}$.
Let $\rho_d (\mathbf{s}) = \prod_{j=1}^d \rho(s_j)$, where:
\begin{align}\label{kl-rho}
\rho(s_j) \equiv \frac{\int K(u) K(u-s_j) du}{\int K^2(u) du },
\end{align}
for each $j$.
We make the following assumptions.
\begin{assum}\label{main-assumption}
Let $d < 4$.
\begin{enumerate}
\item $\mathcal{I} \equiv \prod_{j=1}^d [a_j,b_j]$, where $a_j < b_j, j=1,\ldots,d$, and $\mathcal{I}$ is a strict subset of the support
of $\bm{X}$.
\item The distribution of $\mathbf{X}$ has a bounded Lebesgue density $f_{\mathbf{X}}(\cdot)$ on $\mathbb{R}^d$.
Furthermore, $f_{\mathbf{X}}(\cdot)$ is bounded below from zero with continuous derivatives on $\mathcal{I}$.
\item The density of $U$ is bounded, $\mathbb{E}[U^2|\mathbf{X}=\mathbf{x}]$ is continuous on $\mathcal{I}$, and $\sup_{\mathbf{x} \in \mathbb{R}^d} \mathbb{E}[U^4|\mathbf{X}=\mathbf{x}] < \infty$.
\item $g(\cdot)$ is twice continuously differentiable on $\mathcal{I}$.
\item $\mathbf{K}(\mathbf{s}) = \prod_{j=1}^d K(s_j)$, where
$K(\cdot)$ is a kernel function on $\mathbb{R}$ that has finite support on $[-1,1]$, $\int_{-1}^1 u K(u) du=0$, $\int_{-1}^1 K(u) du=1$,
symmetric around zero, and six times differentiable.
\item $h_n = C n^{-\eta}$, where $C$ and $\eta$ are positive constants such that $\eta < 1/(2d)$ and $\eta > 1/(d+4)$.
\item $\inf_{n \geq 1} \inf_{\mathbf{x} \in \mathcal{I}} s_n(\mathbf{x}) > 0$ and $s_n(\mathbf{x})$ is continuous for each $n \geq 1$.
Furthermore, $\mathbf{x} \mapsto \mathbb{E}\left[ U^2 | \mathbf{X} =\mathbf{x} \right]f_{\mathbf{X}}(\mathbf{x})$ is Lipschitz continuous.
\item There exists an estimator $\hat s^2 (\mathbf{x})$ such that
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{I}} \left| \hat s^2 (\mathbf{x}) - s_n^2(\mathbf{x}) \right| = O_p( n^{-c})
\end{align*}
for some constant $c > 0$.
\item $\max\left\{ (n h_n^d)^{1/2}|\psi (\mathbf{W}_{i},\hat \theta) - \psi (\mathbf{W}_{i},\theta_0) |: i=1,\dots,n \right\} = O_p( n^{-c})$ for some constant $c > 0$.
\end{enumerate}
\end{assum}
Most of the assumptions are standard.
Condition (2) of Assumption \ref{main-assumption} rules out discrete covariates. If all
regressors are discrete, then the estimation problem reduces to a parametric estimation problem.
In this case, one may consider a multiple testing approach as in \cite{Lee:Shaikh:2014} by defining subpopulations with observed cells of discrete covariates. If some covariates are discrete and others are continuous, then one may use a smoothing approach proposed in \cite{Li:Racine:2004}.
Condition (5) assumes that the kernel function has finite support. This assumption is for
the simplicity of the paper and can be dropped at the expense of complicated proofs.
It also assumes that the kernel function is differentiable.
This assumption is crucial and excludes, for example, the uniform kernel.
One of the bandwidth conditions in $h_n$ (that is, $\eta > 1/(d+4)$ in condition (6)) imposes undersmoothing, so that we can ignore the bias asymptotically.
The rule-of-thumb bandwidth proposed in Section \ref{sec:simple:case} satisfies the required rate conditions.
\begin{rem}
Note that $d < 4$ is necessary to ensure that $\eta < 1/(2d)$ and $\eta > 1/(d+4)$ can hold jointly.
It is possible to extend our asymptotic theory to the case that $d \geq 4$ using a higher-order local polynomial estimator under stronger smoothness conditions.
In this paper, we limit our attention to the local linear estimator since we are mainly interested in low dimensional $\bm{x}$.
\end{rem}
\begin{rem}
An estimator of $\hat s^2 (\mathbf{x})$ is readily available. For example, we may consider
\begin{align}\label{sigma-x-est-dim-d}
\hat s^2(\mathbf{x})
=
\frac{\hat{\sigma}^2(\mathbf{x})}{ \hat f_\mathbf{X}(\mathbf{x})}
\int \mathbf{K}^2 \left( \mathbf{u} \right) d \mathbf{u},
\end{align}
where $\hat{f}_\mathbf{X}(\cdot)$ is the kernel density estimator
and
$\hat{\sigma}^2(\mathbf{x})$ is a nonparametric estimator of
${\sigma}^2(\mathbf{x})$ using
$\{ (\hat{U}_i^2, \mathbf{X}_i): i=1,\ldots,n \}$ with
$\hat{U}_i \equiv \psi (\mathbf{W}_i,\hat{\theta} ) - \hat{g}(\mathbf{X}_i)$.
Recall \eqref{sigma-x-est-dim1} for its concrete form for the one-dimensional case.
Alternatively, we may set
\begin{align*}
\hat s^2(\mathbf{x})
= \frac{1}{nh_n^{d}} \sum_{i=1}^n \frac{\hat{U}_i^2}{\hat{f}^2_\mathbf{X}(\mathbf{x})} \mathbf{K}^2 \left( \frac{\mathbf{X} - \mathbf{x}}{h_n} \right).
\end{align*}
For either estimator, it is straightforward to verify condition (8) of Assumption \ref{main-assumption} using the standard results in kernel estimation.
\end{rem}
\begin{rem}
Note that Condition (9) of Assumption \ref{main-assumption} is merely a sufficient (but not necessary) condition.
This condition is satisfied, for example, if
$\| \hat \theta - \theta_0 \| = O_p( n^{-1/2})$,
functions $\beta \mapsto \pi (\mathbf{Z},\beta)$
and
$\alpha_j \mapsto \mu _{j}(\mathbf{Z},\alpha _{j}), j =0,1,$ are
Lipschitz continuous,
$\pi (\mathbf{Z},\beta_0 )$ is bounded between $\epsilon$ and $1-\epsilon$ for some
constant $\epsilon > 0$, provided that
some weak moment conditions on $(Y, \mathbf{Z})$ hold.
\end{rem}
Let $a_n \equiv a_n(\mathcal{I})$ be the largest solution to the following equation:
\begin{align} \label{b1}
\text{mes}(\mathcal{I}) {h_n}^{-d} \lambda^{d/2} (2\pi)^{-(d+1)/2 }
a_n^{d-1} \exp(-a_n^2/2) = 1,
\end{align}
where $\text{mes}(\mathcal{I})$ is the Lebesgue measure of $\mathcal{I}$; that is,
$\text{mes}(\mathcal{I}) = \prod_{j=1}^d (b_j-a_j)$ and:
\begin{equation}\label{def-lambda}
\lambda = \frac{ - \int K(u) K^{\prime \prime }(u) du }{\int K^2(u) du }.
\end{equation}
The following is the main theoretical result of our paper.
\begin{thm}\label{main-uniform}
Let Assumption \ref{main-assumption} hold. Then there exists $\kappa > 0$ such that,
uniformly in $t$, on any finite interval:
\begin{align}\label{main-approx}
\begin{split}
&\mathbb{P}\left( a_n
\left[ \max_{\mathbf{x} \in \mathcal{I}} \left| \frac{\hat{g}(\mathbf{x}) - g(\mathbf{x})}{\hat{s}(\mathbf{x})} \right| - a_n \right] < t \right) = \\
&\exp \left( - 2e^{-t - t^2/2a_n^2} \right)
\sum_{m=0}^{\lfloor (d-1)/2 \rfloor} h_{m,d-1} a_n^{-2m}
\left( 1 + \frac{t}{a_n^2} \right)^{d-2m-1}
+ O( n^{-\kappa}),
\end{split}
\end{align}
as $n \rightarrow \infty$, where $h_{m,d-1} \equiv \frac{(-1)^m (d-1)!}{ m! 2^m (d-2m-1)!}$
and $\lfloor \cdot \rfloor$ is the integer part of a number.
\end{thm}
Notice that the approximation error is of a polynomial rate.
As a result, a critical value based on the leading term of the right-hand side of
\eqref{main-approx} provides a better approximation than one based on the Gumbel approximation.
The result in Theorem \ref{main-uniform} may be of independent interest for constructing the uniform confidence band in nonparametric regression beyond the scope of estimating the CATEF in our context.
\begin{rem}
In a setting different from here, \cite{Lee/Linton/Whang:06} propose analytic critical values based on \cite{Piterbarg:96} in order to test for stochastic monotonicity, compare its performance with the bootstrap critical values in their Monte Carlo experiments, and find that both perform well in finite samples.
However, the discussions in \cite{Lee/Linton/Whang:06} are informal and rely on the results of Monte Carlo experiments without the formal proof of establishing the polynomial approximation error.
\end{rem}
The conclusion of the theorem can be simplified for special cases.
In particular, if $d=1$, then:
\begin{align*}
\mathbb{P}\left( a_n
\left[ \max_{\mathbf{x} \in \mathcal{I}} \left| \frac{\hat{g}(\mathbf{x}) - g(\mathbf{x})}{\hat{s}(\mathbf{x})} \right| - a_n \right] < t \right) = \exp \left( - 2e^{-t - t^2/2a_n^2} \right)
+ O( n^{-\kappa}),
\end{align*}
where $a_n$ is the largest solution to
$\text{mes}(\mathcal{I}) {h_n}^{-1} \lambda^{1/2} (2\pi)^{-1 }\exp(-a_n^2/2) = 1$.
Also, if $d=2$, then:
\begin{align*}
\mathbb{P}\left( a_n
\left[ \max_{\mathbf{x} \in \mathcal{I}} \left| \frac{\hat{g}(\mathbf{x}) - g(\mathbf{x})}{\hat{s}(\mathbf{x})} \right| - a_n \right] < t \right) = \exp \left( - 2e^{-t - t^2/2a_n^2} \right)\left( 1 + \frac{t}{a_n^2} \right)
+ O( n^{-\kappa}),
\end{align*}
where $a_n$ is the largest solution to
$\text{mes}(\mathcal{I}) {h_n}^{-2} \lambda^{2} (2\pi)^{-3/2 }
a_n \exp(-a_n^2/2) = 1$.
\begin{rem}
It is standard to obtain pointwise confidence intervals based on normal approximations.
Recall that our two-sided symmetric uniform confidence band has the form:
\begin{align}\label{pointwise-ci}
\hat{g}(\mathbf{x}) - c(1-\alpha) \frac{\hat{s}(\mathbf{x})}{\sqrt{nh_n^d}}
\leq
g(x)
\leq
\hat{g}(\mathbf{x}) + c(1-\alpha) \frac{\hat{s}(\mathbf{x})}{\sqrt{nh_n^d}},
\end{align}
where $c(1-\alpha)$ is obtained from Theorem \ref{main-uniform}.
To obtain two-sided symmetric pointwise confidence intervals, we just need to replace
$c(1-\alpha)$ with the usual normal critical value $\Phi^{-1}(1-\alpha/2)$, where $\Phi(\cdot)$ is the standard normal cumulative distribution function.
The pointwise confidence interval given in \eqref{pointwise-ci} is different from the one resulting from
\citet[Theorem 2]{Abrevaya-et-al-JBES} in terms of the formula for $\hat{s}(\mathbf{x})$ in \eqref{sigma-x-est-dim-d}. In their case, they need to estimate
$\hat{\sigma}^2(\mathbf{x})$ using
$\{ (\tilde{U}_i^2, \mathbf{X}_i): i=1,\ldots,n \}$ with
\begin{align*}
\tilde{U}_i \equiv \frac{D_i Y_i}{\pi (\mathbf{Z}_i,\hat{\beta} )}- \frac{\left( 1-D_i\right) Y_i}{1-\pi (\mathbf{Z}_i,\hat{\beta})}
\hat{g}(\mathbf{X}_i).
\end{align*}
\end{rem}
\begin{rem}
A one-sided version of the uniform confidence band is readily available. Combining Theorems 14.1 and 14.2 of \cite{Piterbarg:96}
with arguments identical to those used in the proof of Theorem \ref{main-approx-one-sided} yields the following proposition. Under Assumption \ref{main-assumption}, there exists $\kappa > 0$ such that,
uniformly in $t$, on any finite interval:
\begin{align}\label{main-approx-one-sided}
\begin{split}
&\mathbb{P}\left( a_n
\left[ \max_{\mathbf{x} \in \mathcal{I}} \frac{\hat{g}(\mathbf{x}) - g(\mathbf{x})}{\hat{s}(\mathbf{x})} - a_n \right] < t \right) = \\
&\exp \left( - e^{-t - t^2/2a_n^2} \right)
\sum_{m=0}^{\lfloor (d-1)/2 \rfloor} h_{m,d-1} a_n^{-2m}
\left( 1 + \frac{t}{a_n^2} \right)^{d-2m-1}
+ O( n^{-\kappa}),
\end{split}
\end{align}
as $n \rightarrow \infty$.
Note that the only differences between \eqref{main-approx} and \eqref{main-approx-one-sided} are that
(i) there is no absolute value on the left side of the equation in \eqref{main-approx-one-sided}
and (ii) there is no factor 2 inside the exponential function in \eqref{main-approx-one-sided}.
Hence, for example, if $d=1$, then:
\begin{align*}
\mathbb{P}\left( a_n
\left[ \max_{\mathbf{x} \in \mathcal{I}} \frac{\hat{g}(\mathbf{x}) - g(\mathbf{x})}{\hat{s}(\mathbf{x})} - a_n \right] < t \right) = \exp \left( - e^{-t - t^2/2a_n^2} \right)
+ O( n^{-\kappa}).
\end{align*}
\end{rem}
\subsection{Construction of critical values}
We use the leading term on the right-hand side of \eqref{main-approx}
as a distribution-like function to construct a uniform confidence band.
For example, if $d=1$, we may construct a critical value $c(1-\alpha)$ that satisfies:
\begin{align*}
F_{n,1}(c) \geq 1 - \alpha,
\end{align*}
where $F_{n,1}(t) \equiv \exp \left( - 2e^{-t - t^2/2a_n^2} \right)$.
This yields the critical value presented in the Algorithm of Section \ref{sec:simple:case}.
Similarly, if $d=2$, we can use:
\begin{align*}
F_{n,2}(c) \geq 1 - \alpha,
\end{align*}
where $F_{n,2}(t) \equiv \exp \left( - 2e^{-t - t^2/2a_n^2} \right)\left( 1 + \frac{t}{a_n^2} \right)$.
In finite samples, it might be useful to impose monotonicity of $F_{n,j}(\cdot)$
by rearrangement
(see, e.g., \cite{CF-VG:09}).
\begin{rem}
Theorem \ref{main-uniform} implies that:
\begin{align*}
&\lim_{n \rightarrow \infty} \mathbb{P}\left( a_n
\left[ \max_{\mathbf{x} \in \mathcal{I}} \left| \frac{\hat{g}(\mathbf{x}) - g(\mathbf{x})}{\hat{s}(\mathbf{x})} \right| - a_n \right] < t \right) = \exp \left( - 2e^{-t } \right).
\end{align*}
Thus, one may construct analytical critical values based on the Gumbel distribution. However, this approximation is accurate only up to the logarithmic rate in view of Theorem \ref{main-uniform}.
\end{rem}
\section{Monte Carlo Experiments}\label{sec:MC}
In this section, we present the results of Monte Carlo experiments.
These experiments are conducted to see the finite sample performances of the proposed doubly robust estimator and the proposed uniform confidence band.
The simulations are conducted by R 3.3.1 with Windows 10.
The number of replications is $5000$.
\subsection{Data generating process}
The data generating process follows the potential outcome framework.
The notations for the variables are the same as those used in the theoretical part of the paper.
We consider cases with $p=10, 30$ and $N= 500, 2000$.
The data generating process is the following.
The vector of covariates $\mathbf{Z}= (X_1, \dots, X_p)^{\top} $ is generated by:
$\mathbf{Z} \sim N(0,I_p)$, where $I_p$ is the $p$-dimensional identity matrix.
The potential outcomes are generated by:
\begin{align*}
Y_1 = 10 + \sum_{k=1}^p \frac{1}{\sqrt{p}} X_k + v, \quad Y_0 =0,
\end{align*}
where $v \sim N(0, 1)$ and $v$ is independent of $\mathbf{Z}$.
The treatment status $D$ is generated by:
\begin{align*}
D =\mathbf{1} \left\{ \Lambda \left(\sum_{k= p/2}^p \frac{1}{\sqrt{p/2}}X_k \right) > U\right\},
\end{align*}
where $U \sim U[0,1]$, $U$ is independent of $(\mathbf{Z}^{\top}, v)$
and $\Lambda$ is the logistic function.
Thus, the propensity score is $\pi(\mathbf{Z}) = \Lambda \left(\sum_{k= p/2}^p X_k/ \sqrt{p/2}\right)$.
The observed outcome is $Y=DY_1$.
The parameter of interest is the CATEF for $\mathbf{X} =X_1$.
In our specification, the CATEF
can be written as:
\begin{align*}
CATEF (x_1) = 10 + x_1/\sqrt{p} .
\end{align*}
We examine the performance of various statistical procedures regarding this CATEF.
\subsection{Model specification}
To estimate and conduct statistical inferences on $CATEF(x_1)$ using our doubly robust procedure, we need to specify a model for the regression $\mu_j (\mathbf{z})$ for $j=0,1$ and a model for the propensity score $\pi (\mathbf{z})$.
We consider two regression models and two propensity score models.
One of two models is correctly specified, but the other model is misspecified.
We note that our doubly robust procedure is predicted to work well provided that at least one of the regression model and the propensity score model is correctly specified.
We first discuss the model specifications for the regression part.
The first regression model is:
\begin{align*}
\mu_1 (\mathbf{z}, \alpha_1) = \alpha_{10} + \sum_{k=1}^p \alpha_{1k} X_k, \quad
\mu_0 (\mathbf{z}, \alpha_0) = \alpha_{00} + \sum_{k=1}^p \alpha_{0k} X_k .
\end{align*}
This model is correctly specified.
The coefficients are estimated by OLS using $(1, X_1, \dots, X_p)$ as the explanatory variable.
The second regression model, which is misspecified, is:
\begin{align*}
\mu_1 (\mathbf{z}, \alpha_1) = \alpha_{10} + \sum_{k=1}^{p/2} \alpha_{1k} X_k ,\quad
\mu_0 (\mathbf{z}, \alpha_0) = \alpha_{00} + \sum_{k=1}^{p/2} \alpha_{0k} X_k .
\end{align*}
The model is estimated by OLS using $(1, X_1, \dots, X_{p/2})$ as explanatory variables.
This model is misspecified because it suffers from sample selection bias introduced by omitting the second half of the regressors which affects the treatment status.
We also consider two models for propensity score.
The model for propensity score is:
\begin{align*}
\pi (\mathbf{z}, \beta) = \Lambda \left(\beta_0 + \sum_{k=1}^p \beta_k X_k\right).
\end{align*}
The misspecified model is:
\begin{align*}
\pi(\mathbf{z},\beta) = \Lambda \left(\beta_0 + \sum_{k=1}^{p/2} \beta_k X_k \right).
\end{align*}
Similarly to the case of the regression part, misspecification is introduced by omitting the second half of the regressors.
The models for propensity score are estimated by maximum likelihood.
We estimate $CATEF(x_1)$ for $x_1 \in \{ -1, -0.5, 0, 0.5, 1\}$ and compute the mean bias (``MEAN''), standard deviation (``SD''), the average of standard error for $\widehat{CATEF} (x_1)$ (``ASE''), and the root mean squared error (``RMSE'').
The local linear regression is conducted with the Gaussian kernel, and the preliminary bandwidth ($\hat h$ in Algorithm (1)) is chosen by the method of \cite{Ruppert:Sheather:Wand:1995}.
We also compute the ``BIAS'', ``SE'' and ``RMSE'' of the corresponding inverse probability weighting estimators and the regression adjustment estimators.
Note that the difference between the proposed method and those alternative methods arises only in the estimation of $\psi (\mathbf{W}, \theta_0)$ and the other steps are the same.
We examine the coverage probability of the uniform confidence band for $CATEF (x_1)$ for the range $-1 \le x_1 \le 1$.
The nominal coverage probabilities that we consider are 99\%, 95\% and 90\%.
We compute the empirical coverage (``CP''), the mean critical value (``Mcri''), and the standard deviation of critical value (``Sdcri'').
We also compute the coverage probabilities of the confidence band based on the critical values computed by the Gumbel approximation (``GCP'').
\subsection{Results}
Tables \ref{tab-est-k10} and \ref{tab-est-k30} summarize the results on the properties of the estimators.
In both tables, DR refers to our doubly robust method, whereas IPW and RA correspond to the inverse probability weighting and regression adjustment methods, respectively.
The proposed doubly robust estimator of the CATEF works well in finite samples.
As the theory indicates, the proposed estimator exhibits small bias provided that at least one of the regression model and the propensity score model is correctly specified.
We find that the regression adjustment estimator is very precise when the regression model is correctly specified.
However, it suffers from substantial bias when the regression model is misspecified.
The inverse probability weighting estimator also suffers from model misspecification.
Moreover, its standard deviation is much larger than those of the doubly robust and regression adjustment estimators.
When both models are misspecified, all three estimators suffer from heavy bias.
The inverse probability weighting estimator has the largest RMSE because its distribution is more diverse than those of the other two estimators.
All the estimators have larger standard deviations when $x=1$ and $x=-1$ compared to those in other points.
This is because the number of observations around $x=1$ or $x=-1$ is expected to be smaller than that around, for example, $x=0$ which is the center of the distribution.
On the other hand, the magnitude of the bias does not vary much across data points.
The standard error for the proposed doubly robust estimator is slightly smaller than the standard deviation, but the difference is not large.
Tables \ref{tab-cb-k10} and \ref{tab-cb-k30} summarize the finite sample properties of the proposed
uniform confidence band.
The results show that our uniform confidence band has a reasonably good coverage property provided that one of the models is correctly specified.
When both models are misspecified, the size distortion is heavy.
We observe that the size distortion is heavier when the regression model is misspecified than that in the case of propensity score misspecification.
This result indicates that we should carefully model the regression part in order to obtain reliable confidence bands.
The average values of the 95\% critical values are around 2.75. Because the pointwise critical value is 1.96 and is much smaller than the uniformly valid critical value, it demonstrates the importance of the uniform property of confidence band.
The standard deviations of the critical values are small because they change only if the bandwidth changes.
The confidence band based on the Gumbel approximation is very conservative.
The results of the Monte Carlo simulation confirm that the proposed doubly robust estimator indeed works well in finite samples provided that one of the regression and propensity score models is correctly specified.
The proposed uniform confidence band also has good coverage properties.
\section{An Empirical Application}\label{sec:EA}
We apply our uniformly valid confidence band for the CATEF for the effect of maternal smoking
on birth weight where the argument of the CATEF is the mother's age.
Our aim here is to illustrate our confidence band in comparison with alternative confidence bands.
We first discuss the background of this application and the datasets used.
We use two different data sets: the dataset from Pennsylvania and that from North Carolina.
We then compute various confidence bands for the CATEF and discuss the results.
While the purpose of this application is to illustrate our uniformly valid confidence band and not to present new insights on the effect of smoking, it is still informative to discuss the background of this application.
Many studies document that low birth weight is associated with prolonged negative effects on health and educational or labor market outcomes throughout life, although there has been a debate over its magnitude. See, e.g., \cite{AlmondCurrie11} for a review.
Maternal smoking is considered to be the most important preventable negative cause of low birth weight \citep{Kramer87}.
There are many studies that evaluate the effect of maternal smoking on low birth weight \citep{AlmondCurrie11}.
The program evaluation approach is employed by, for example,
\cite{AlmondChayLee05}, \cite{daVeigaWilder08} and
\cite{WalkerTekinWallace09}, and panel data analysis is carried out by
\cite{Abrevaya06} and \cite{AbrevayaDahl08}.
Here, we are interested in how the effect of smoking changes across different age groups of mothers.
\cite{WalkerTekinWallace09} examine whether the effect of smoking is larger for teen mothers than for adult mothers and find mixed evidence.
\cite{Abrevaya-et-al-JBES} also consider this problem in their application.
\subsection{Pennsylvania data}
The first dataset consists of observations from white mothers in Pennsylvania in the USA.
The dataset is an excerpt from \cite{Cattaneo10} and is obtained from the STATA website (``\url{http://www.stata-press.com/data/r13/cattaneo2.dta}'').
Note that the dataset was originally used in \cite{AlmondChayLee05}.
We restrict our sample to white and non-Hispanic mothers, and the sample size is 3754.
The outcome of interest ($Y$) is infant birth weight measured in grams.
The treatment variable ($D$) is a binary variable that is equal to 1 if the mother smokes and 0 otherwise.
The set of covariates $\mathbf{Z}$ includes the mother's age, an indicator variable for alcohol consumption during pregnancy, an indicator for the first baby, the mother's educational attainment, an indicator for the first prenatal visit in the first trimester, the number of prenatal care visits, and an indicator for whether there was a previous birth where the newborn died.
We are interested in how the effect of smoking varies across different values of the mother's age.
Therefore, $\mathbf{X}$ is mother's age in this application.
To estimate the CATEF, we use linear regression models for the regression part and a logit model for propensity score.
The explanatory variables used in the regression models and the logit model
consist of all the elements of $\mathbf{Z}$, the square of the mother's age, and the interaction terms between the mother's age and all other elements of $\mathbf{Z}$.
We estimate the CATEF in the interval between ages 15 and 35.
We compute the following three 95\% confidence bands for the CATEF.
``Our CB'' is the confidence band proposed in this paper.
Because $\mathbf{X}$ is univariate in this application, we follow the algorithm in Section 3.
We use the Gaussian kernel.
The preliminary bandwidth ($\hat h$) is chosen by the method of \cite{Ruppert:Sheather:Wand:1995}.
``Gumbel CB'' is the confidence band in which $c(1-\alpha)$ in the algorithm is replaced by that based on the Gumbel approximation (see Remark \ref{rem-gumbel}).
``PW CB'' is a pointwise valid confidence band where we replace $c(1-\alpha)$ in the algorithm by the corresponding value from the standard normal distribution (i.e., $1.96$).
This provides a valid confidence interval for each point of the CATEF.
However, its uniform coverage rate would be smaller than 95\%.
Figure 1 plots the estimated CATEF and the three 95\% confidence bands for the range between
15 and 35 years of age.
The figure also contains the average treatment effect estimate (AIPW estimate) for a reference.
The widths of the three confidence bands are substantially different.
The confidence band based on the Gumbel approximation provides the widest band and may not be very informative.
The confidence band that is valid only in a pointwise sense gives the narrowest band.
This band is not uniformly valid and so may provide misleading information about the CATEF. On the other hand, this provides valuable information if we are interested at a particular point of the CATEF.
The confidence band we propose lay between ``Gumbel CB'' and ``PW CB''.
While this band is wider than ``PW CB'', it is much narrower than ``Gumbel CB'' and is uniformly asymptotically valid.
We see from this figure that our confidence band is informative while being uniformly valid.
The estimated CATEF is decreasing from 15 to around 25 years of age. It is rather stable for the range above 25 years of age.
All confidence bands indicate that the CATEF is estimated imprecisely near the ends of the range.
Nonetheless, the estimated CATEF indicates that smoking may not have a strong impact when the mother is young.
The CATEF is estimated relatively precisely in the middle of the range.
For the range between 20 and 30 years of age, even the band based on the Gumbel approximation, which is the widest, does not contain 0.
This result provides robust evidence that smoking has a negative impact on birth weight at least for mothers who are 20 to 30 years old.
In this particular dataset, the statistical evidence against a constant smoking effect is somewhat weak.
The confidence band that is valid only in a pointwise sense may provide an impression that the smoking effect depends on the mother's age.
However, the uniformly valid confidence band that we propose marginally contains the straight line
that is equal to the ATE estimate.
This result illustrates that there is a caveat when we use pointwise confidence intervals, as well as the importance of using uniformly valid confidence bands.
\subsection{North Carolina data}
The second dataset is based on the records between 1988 and 2002 by the North Carolina State Center Health Services.
This dataset is used in \cite{Abrevaya-et-al-JBES} and obtained from Robert Lieli's website (``\url{http://www.personal.ceu.hu/staff/Robert_Lieli/cate-birthdata.zip}).
We restrict our sample to white and first--time mothers, and the sample size is 433,558.
As in the case of the Pennsylvania data, the outcome is infant birth weight measured in grams
and the treatment variable is an indicator for smoking status.
The set of covariates $\mathbf{Z}$ includes those used in the analysis of the Pennsylvania data, except an indicator for the first baby because we focus on first--time mothers, and in addition, it includes indicators for gestational diabetes, hypertension, amniocentesis and ultra sound exams.
Again, $\mathbf{X}$ is mother's age in this application.
The specification for the estimation of the CATEF is the same as before.
The purpose of using this much larger dataset is to examine the effect of the sample size.
Our method involves nonparametric kernel regression and it might require a large sample size to yield a reliable result.
For example, the result from the Pennsylvania data indicates that the effect of smoking is very small for very young mothers.
One might argue that such a result is an artifact of small sample size.
The other issue is that the confidence bands obtained using the Pennsylvania data are somewhat wide.
We hope that using this larger dataset provides us with narrower confidence bands and more informative statistical results.
Figure 2 plots the estimated CATEF and the three 95\% confidence bands for the range between
15 and 35 years of age.
Note that the scale of the vertical axis is different from Figure 1.
We now obtain much narrower confidence bands.
The widths of the three (uniform, point-wise and Gumbel) confidence bands
are still different.
The estimated CATEF for young mothers is negative and statistically different from 0.
The previous result that it is close to 0 may be considered as an artifact of small sample size.
The estimated CATEF is decreasing from around 17 to around 29 years of age.
For the range above 30 years of age, we obtain relatively wide confidence bands.
We reject the null hypothesis of no effect of smoking on birth weights uniformly over 15-35 years of age.
These confidence bands do not support the hypothesis that the CATEF is constant because the ATE line exceeds the confidence bands.
One might argue that the difference in the results may stem from the fact that the North Carolina data contains richer information and we use a larger set of covariates.
We reexamine the North Carolina data based on the same set of covariates as that for the Pennsylvania data, except an indicator for the first baby.
Figure 3 plots the estimated CATEF and confidence bands obtained using this set of covariates.
The results in Figure 3 are qualitatively very similar to those in Figure 2.
We thus believe that the difference between the results from the Pennsylvania data and the North Carolina data are not from the difference in the covariates but from the difference in the sizes of these two samples.
We thus interpret our findings to indicate that the different results come from the difference in sample size yet our confidence bands reasonably quantify the uncertainty from small sample size.
While two data-sets yield different estimates of CATEF, the confidence bands from the Pennsylvania data include the estimated CATEF and the confidence bands from the North Carolina data.
While we use the same data set as that used in \cite{Abrevaya-et-al-JBES}, it is somewhat difficult to compare their results with ours because of differences in the implementations.
In particular, the bandwidths are very different.
Our choice of bandwidth is around 0.2, while theirs are between 1.4--11.2.
Nonetheless, we make several remarks.
Using small bandwidths (1.4 and 2.8), \cite{Abrevaya-et-al-JBES} observe almost no effect for young mothers and a large negative effect for 25--30 years old mothers.
We do not observe such a large difference in the effect across different age groups.
Our confidence band is as tight as their confidence band obtained with bandwidth 11.2 even though we use a much smaller bandwidth and our confidence band is uniform.
This is possibly because we use an AIPW method which yields a more efficient estimate than an IPW method does.
\section{Conclusion}\label{sec:C}
In this paper, we propose a doubly robust method for estimating the CATEF.
We consider the situation where a high-dimensional vector of covariates is needed for identifying the average treatment effect but the covariates of interest are of much lower dimension.
Our proposed estimator is doubly robust and does not suffer from the curse of dimensionality.
We propose a uniform confidence band that is easy to compute, and we illustrate its usefulness via Monte Carlo experiments and an application to the effects of smoking on birth weights.
There are a few topics to be explored in the future.
First, it would be useful to consider the issue of asymptotic biases of the proposed estimator without relying on undersmoothing. For example, it might be possible to extend the approach
of \cite{hall2013} that avoids undersmoothing for our purposes.
Second, it would be an interesting exercise to develop a method for estimating the quantile treatment effects conditional on covariates.
Third, it is possible to extend our approach to the local average treatment effect.
As mentioned in the Introduction,
\cite{Ogburn-et-al-JRSSB} consider conditioning on $\mathbf{Z}$ to achieve identification, but they estimate the local average treatment effect, say LATE($\mathbf{X}$), as a function of $\mathbf{X}$. However, their specification of LATE($\mathbf{X}$) is parametric.
Our approach can be adapted to specify LATE($\mathbf{X}$) nonparametrically and to develop a corresponding uniform confidence band.
Fourth, this paper does not cover marginal treatment effects that can be identified using the method of local instrumental variables developed by \citet{Heckman:Vytlacil:99,Heckman:Vytlacil:05}. It would be interesting to develop a uniform confidence band for the marginal treatment effects.