\begin{document}

\def\spacingset#1{\renewcommand{\baselinestretch}
{#1}\small\normalsize} \spacingset{1}

\if00
{
  \title{ 
 \bf Robust causal inference with \\ continuous instruments using the \\ local instrumental variable curve}
  \author{Edward H. Kennedy
    \thanks{Edward Kennedy is Assistant Professor of Statistics, Carnegie Mellon University, Pittsburgh, PA 15217 (e-mail: edward@stat.cmu.edu). Scott Lorch is Associate Professor of Pediatrics at the University of Pennsylvania School of Medicine and Attending Neonatologist at The Children's Hospital of Philadelphia. Dylan Small is Professor of Statistics, University of Pennsylvania. The authors gratefully acknowledge support from NIH grant R01-DK090385 (Kennedy) and NSF grant SES-1260782 (Small).}\hspace{.2cm}
    \\   \date{}
    Department of Statistics, Carnegie Mellon University \\ \\
    Scott A. Lorch \\
    Department of Pediatrics, Perelman School of Medicine, \\ University of Pennsylvania \\ \\
    Dylan S. Small \\
    Department of Statistics, The Wharton School, University of Pennsylvania }
  \maketitle
  \setcounter{page}{0}
  \thispagestyle{empty}

} \fi

\if10
{
  

  

  

  \begin{center}
    {\LARGE\bf Robust causal inference with \\ 
 continuous  instruments  using the \\ 
  local instrumental variable curve}

  \end{center}
  \setcounter{page}{0}
  

} \fi

\begin{abstract}
Instrumental variables are commonly used to estimate effects of a treatment afflicted by unmeasured confounding, and in practice instruments are often continuous (e.g., measures of distance, or treatment preference). However, available methods for continuous instruments have important limitations: they either require restrictive parametric assumptions for identification, or else rely on modeling both the outcome and treatment process well (and require modeling effect modification by all adjustment covariates). \textcolor{black}{In this work we develop the first semiparametric doubly robust  estimators of the local instrumental variable effect curve}, i.e., the effect among those who would take treatment for instrument values above some threshold and not below. In addition to being robust to misspecification of either the instrument or treatment/outcome processes, our approach also incorporates information about the instrument mechanism and allows for flexible data-adaptive estimation of effect modification. We discuss asymptotic properties under weak conditions, and use the methods to study infant mortality effects of neonatal intensive care units with high versus low technical capacity, using travel time as an instrument. 
\end{abstract}

\noindent
{\it Keywords:} causal inference, complier average treatment effect, cross-validation, doubly robust, semiparametric theory. \\

\thispagestyle{empty}

\spacingset{1.45} 
\spacingset{1}

\section{Introduction}
\label{sec:intro}

Instrumental variables (IVs) are commonly used to estimate effects of treatments that are afflicted by unmeasured confounding. Instruments are special variables that influence treatment, but are themselves unconfounded and do not directly affect outcomes, allowing the recovery of some causal information from data that might otherwise be unusable. In practice, instruments are often continuous (e.g., measures of distance, or treatment preference), but most available methods only consider instruments that are discrete (and typically binary). Further, methods that do allow for continuous IVs have important limitations.

Classical IV methods (e.g., standard two-stage least squares), which were developed in a structural equation model framework, allow for continuous instruments but require strong parametric assumptions for identification, assume that treatment effects do not vary across units, and also require correct parametric models for at least how the outcome process depends on covariates and instruments \autocite{wooldridge2010econometric, okui2012doubly}. Alternatively, Robins and others  \autocite{robins1989analysis, robins1994correcting, hernan2006instruments, tan2010marginal, okui2012doubly} developed approaches in the potential outcomes framework that can also handle continuous instruments, but which allow heterogeneous treatment effects, and also permit doubly robust covariate adjustment \autocite{robins2001inference,van2003unified,bang2005doubly}. Doubly robust IV methods are consistent as long as either the instrument mechanism or the treatment/outcome mechanisms are correctly modeled (not necessarily both), and they can also yield fast root-n convergence rates and inference even when using flexible nonparametric methods for covariate adjustment. However, the methods developed in this framework still require parametric assumptions for identification; they typically target treatment effects on the treated, and achieve identification with dimension-reducing parametric assumptions that restrict how heterogeneous treatment effects can be. As noted for example by \textcite{tchetgen2013alternative}, this kind of approach is problematic because a priori information about the parametric form of underlying causal structure is rarely available, and misspecification could lead to large biases that cannot be detected with data.

An alternative approach is to replace dimension-reducing homogeneity assumptions with a monotonicity assumption  \autocite{robins1989analysis, imbens1994identification}, which rules out the possibility that any units would respond oppositely to encouragement from the instrument. In other words, there can be units who are encouraged by the instrument, as well as units who do not respond at all to the instrument, but there cannot be units who defy encouragement from the instrument. For example, in the binary instrument case, there can be units who take treatment if and only if they receive the instrument, as well as units who always or never take treatment; however, there cannot be units who take control if the instrument is received but take treatment if not.

This assumption is often plausible in practice, and also permits nonparametric identification of causal effects among compliers (i.e., those who do respond to encouragement from the instrument). However, monotonicity is usually framed in terms of binary instruments \autocite{imbens1994identification, abadie2003semiparametric, tan2006regression, ogburn2015doubly}. An important exception is a strand of work that has focused on estimating local IV (LIV) curves, i.e., effects among units who would comply right at a given threshold value of the instrument \autocite{heckman1997instrumental, heckman1999local,glickman2000derivation, heckman2005structural}.

This literature on LIV approaches arose out of a latent index or selection model framework \autocite{vytlacil2002independence}, and is unique in allowing for continuous instruments while still permitting nonparametric identification. However there are important limitations. {First, available approaches for estimating the LIV curve rely on modeling how both the treatment and outcome depend on covariates and instrument \autocite{basu2007use, carneiro2010estimating}, and typically use restrictive parametric models. This is problematic since parametric models are often relied upon based on convenience, rather than real substantive knowledge, and can yield severe bias if misspecified. Conversely, fully nonparametric approaches are sensitive to the curse of dimensionality and typically yield estimators with slow rates of convergence, as well as little hope for centered confidence intervals without impractical undersmoothing {(we refer to Section 5.7 of \textcite{wasserman2006all} for details)}. Further, in the IV setting there may be some information available about how the instrument depends on covariates (e.g., about the density of the instrument given covariates), but this is not incorporated in approaches that rely solely on treatment and outcome models, whether parametric or nonparametric. Our paper solves these problems with a semiparametric doubly robust approach that can attain parametric rates of convergence, even while allowing flexible nonparametric estimation of nuisance functions.}

Second, available LIV estimands are fully conditional on all measured covariates, even though in many cases effect modification is not of particular scientific interest, or else it is only of interest for a small subset of covariates. Marginal effects {can often be estimated more robustly at faster rates of convergence,} and are often more closely tied to scientific questions. \textcite{van2003unified} point out that using fully conditional effects puts us at the whim of whatever confounders happen to arise in the dataset at hand, whereas marginal effects allow framing scientific questions a priori. Marginalizing currently available fully conditional estimators leads to awkward and uninterpretable models, as discussed in the local average treatment effect (LATE) setting by \textcite{ogburn2015doubly}. {However, in contrast to the LATE setting, currently available methods cannot be adapted for doubly robust estimation of the LIV curve, even in the simpler case where effect modifiers are discrete or not present. Our paper solves these problems by developing methods for working models of a marginal version of the LIV curve itself, allowing for arbitrary effect modification.  Our direct modeling approach eases interpretability by allowing analysts to incorporate background knowledge on the actual parameter of interest, without having to specify models for nuisance quantities that are not of direct scientific interest.}

In addition to the above, our work makes several other important advances. {Importantly, we use empirical process theory and sample splitting to derive asymptotic properties of our approach under weak conditions, which allow for flexible data-adaptive estimation of nuisance functions in the presence of complex high-dimensional confounding. We also develop a doubly robust cross-validation approach for model selection in high-dimensional settings, which is crucial for learning the LIV curve from data.} Finally, we explore finite-sample properties via simulation, and implement our methods to study effects of high-level neonatal intensive care units (NICUs) on infant mortality, using travel time as an instrument.

\section{Preliminaries}

\subsection{Data \& Notation}
\label{4data}

Suppose we observe an independent and identically distributed sample $(\mathbf{O}_1,...,\mathbf{O}_n)$ with $\mathbf{O}=(\mathbf{X},Z,A,Y)$, where $\mathbf{X}$ is a vector of covariates, $Z$ is a continuous instrument for a binary treatment $A$, and $Y$ is some real-valued outcome of interest. The covariates $\mathbf{X}=(\mathbf{V},\mathbf{W})$ are partitioned into potential effect modifiers of interest $\mathbf{V}$ and other covariates $\mathbf{W}=\mathbf{X} \setminus \mathbf{V}$ not of interest but for which adjustment is still necessary. The choice of $\mathbf{V}$ is based purely on the scientific question, so that if effect modification is not of interest one can simply select $\mathbf{V}=\emptyset$. We characterize causal effects using potential outcome notation \autocite{rubin1974estimating}, and so let $Y^a$ (and $Y^{za}$) denote the potential outcomes that would have been observed had treatment level $A=a$ (and instrument level $Z=z$) been received. Similarly we let $A^z$ denote the potential treatment that would have been observed under instrument level $Z=z$. A directed acyclic graph showing the data structure is given in Figure~\ref{fig:dag}. 

\begin{figure}
\caption{Directed acyclic graph showing covariates $\mathbf{X}
\label{fig:dag}
\end{figure}

We let $P$ denote the distribution of $\mathbf{O}$, with density with respect to some dominating measure given by $p(\mathbf{o})=p(y \mid \mathbf{x},z,a) p(a \mid \mathbf{x},z) p(z \mid \mathbf{x}) p(\mathbf{x})$. In general we write the density of a variable $T$ under $P$ evaluated at $z$ as $p(T=z)$, except when there is no ambiguity (e.g., $p(t)$ is the density of $T$ at $t$), and we use $\text{supp}(T)$ for the support of $T$. Finally we use some additional notation to simplify the presentation. Specifically we let $\pi(z \mid \mathbf{x})=p(Z=z \mid \mathbf{X}=\mathbf{x})$ denote the density of the instrument given covariates (i.e., instrument propensity score), $\mu(\mathbf{x},z)=\mathbb{E}(Y \mid \mathbf{X}=\mathbf{x}, Z=z)$ and $\lambda(\mathbf{x},z)=\mathbb{E}(A \mid \mathbf{X}=\mathbf{x}, Z=z)$ denote the outcome and treatment regression functions, respectively, with marginalized versions given by $m(z,\mathbf{v})=\mathbb{E}\{ \mu(\mathbf{X},z) \mid \mathbf{V}=\mathbf{v}\}$ and $\ell(z,\mathbf{v})= \mathbb{E}\{ \lambda(\mathbf{X},z) \mid \mathbf{V}=\mathbf{v}\}$. We let $\mathbb{P}_n$ denote the empirical measure so that empirical averages can be written as $n^{-1} \sum_i f(\mathbf{O}_i)=\mathbb{P}_n\{f(\mathbf{O})\}$. The notation $|| \cdot ||$ denotes the Euclidean norm $|| \mathbf{o}ldsymbol\beta ||=(\mathbf{o}ldsymbol\beta^{ \mathrm{\scriptscriptstyle T} } \mathbf{o}ldsymbol\beta )^{1/2}$, and $||f||_2=\{\int f(\mathbf{o})^2 \ dP(\mathbf{o})\}^{1/2}$ denotes the $L_2(P)$ norm .

\subsection{Monotonicity}

Before defining the causal estimand of interest and considering identifying assumptions, it will be helpful to discuss the concept of monotonicity, which was introduced in various forms by \textcite{robins1989analysis} and \textcite{imbens1994identification}, among others. In a classical binary instrument setting (where supp$(Z)=\{0,1\}$), monotonicity can be stated succinctly as
$$ A^1 \geq A^0 \text{ with probability one}. $$
Monotonicity rules out the possibility that there are troublesome units in the population with $A^0=1$ but $A^1=0$. Such units are called `defiers' since they take treatment $A=1$ when not encouraged by the instrument (i.e., when $Z=0$), but take control $A=0$ when they are in fact encouraged (i.e., when $Z=1$). Thus monotonicity ensures that the population only comprises never-takers ($A^0=A^1=0$), always-takers ($A^0=A^1=1$), and compliers ($A^0=0, A^1=1$). Monotonicity can often be a reasonable assumption in practice, but not always; it has been discussed extensively in previous work, particularly for binary instruments (see \textcite{imbens2014instrumental} and discussion for a nice overview).

A natural way to extend monotonicity to the continuous instrument setting is as follows.

\begin{assumption}[Monotonicity] \label{ass:4mono}
If $z' > z$ then $A^{z'} \geq A^z$ with probability one.
\end{assumption}

Under Assumption \ref{ass:4mono}, no unit would ever change from treatment to control with an increase in the instrument value; increasing the instrument can either encourage treatment over control or have no effect at all, but it cannot discourage treatment relative to lesser instrument values. Thus the population still comprises never-takers, always-takers, and compliers, but with continuous instruments the compliers can be further partitioned into compliers at given instrument values. In particular, a complier at $Z=z$ would be a unit for which $A^z=1$ but $A^{z-\delta}=0$ for any $\delta>0$. {We use a population version of monotonicity for simplicity, but only a covariate-specific conditional version is necessary when targeting conditional treatment effects \autocite{tan2006regression}.} The above continuous version of monotonicity has been employed and discussed  by \textcite{glickman2000derivation} and  \textcite{vytlacil2002independence}, for example. Importantly, these authors showed that (when coupled with standard identifying assumptions to be discussed shortly) the above monotonicity assumption is equivalent to the following latent threshold model. 

\begin{assumptionprime}[Latent Threshold] \label{ass:4thresh}
$A^z = \mathbbold{1}(z \geq T)$ for an unobserved random threshold $T$.
\end{assumptionprime}

Under the latent threshold model, each complier has some instrument value at which they are encouraged to take treatment, while for any lesser value they would take control. Larger values of the threshold $T$ indicate units that are less willing to take treatment, i.e., less susceptible to encouragement by the instrument. We can thus define the latent threshold $T$ as 
$$ T = \begin{cases} -\infty & \text{if $A^z=1$ for all $z$ (always-takers)}  \\ \inf \{ z  : A^z =1 \} & \text{if $A^{z'}>A^z$ for some $z'>z$ (compliers)} \\ \infty & \text{if $A^z=0$ for all $z$ (never-takers).} \end{cases}$$
We refer to \textcite{vytlacil2002independence} for further discussion and detail.

\subsection{Estimand \& Identification}

In this paper our goal is estimation and inference for the local instrumental variable (LIV) curve, which we define as
\begin{equation} \label{eq:liv}
\gamma(t,\mathbf{v}) = \mathbb{E}(Y^1 - Y^0 \mid T=t, \mathbf{V}=\mathbf{v}) .
\end{equation}
This is the average treatment effect among those with latent threshold $T=t$ (and value $\mathbf{V}=\mathbf{v}$ for some arbitrary baseline covariate subset), i.e., the effect among those units with $\mathbf{V}=\mathbf{v}$ who would be encouraged to take treatment right when the instrument passes $Z=t$ but not for lesser values. A fully conditional version of the LIV curve with $\mathbf{V}=\mathbf{X}$ was proposed in the latent index or selection model framework by \textcite{heckman1997instrumental}, and discussed in detail by \textcite{heckman1999local} and \textcite{heckman2005structural}. { Note that the variables $\mathbf{V}$ are arbitrary, and no causal assumptions need to be made about them; this is because we are estimating effect modification (how effects vary with $\mathbf{V}$) and not mediation (how effects would change if we set $\mathbf{V}$ to certain values). }

{We study the LIV curve for continuous IVs since such IVs are very common in practice and since novel theory is required; an analogous LIV effect  $\mathbb{E}(Y^1 - Y^0 \mid A^t > A^{t-1}, \mathbf{V})$  for discrete multivalued instruments could be studied by adapting results from \textcite{tan2006regression, ogburn2015doubly}. However, for continuous instruments, the LIV effect curve is substantially different from the local average treatment effect (LATE) studied by \textcite{angrist1996identification}, \textcite{tan2006regression}, \textcite{ogburn2015doubly} and others. Specifically, the continuous IV extension of the LATE is the three-dimensional effect surface $\mathbb{E}(Y^1-Y^0 \mid A^z>A^{z'}, \mathbf{V}=\mathbf{v})$ for values $z>z'$, which represents the effect among those who would take treatment at $Z=z$ but not $Z=z'$. In contrast, for continuous instruments, the LIV parameter in \eqref{eq:liv} is an arguably more interpretable (and easier to model) curve representing the effect among those ``at the margin'' who would take treatment at $Z=t$ but not $Z<t$. These and other differences between the LATE and LIV are discussed in more detail by \textcite{heckman1999local, heckman2005structural}. Further, even more substantial differences arise after identification, as discussed below. }

\begin{remark}
In Heckman's selection model framework, the LIV curve in \eqref{eq:liv} with $\mathbf{V}=\mathbf{X}$ was termed the ``marginal treatment effect'', and its observed data counterpart the ``local IVs'' estimand (after employing identifying assumptions). We use the latter in both cases, since ``marginal'' is often used to mean ``averaged'' instead of ``at the margin''.
\end{remark}

For identification we rely on standard IV assumptions, which have been employed for example by \textcite{angrist1996identification}, \textcite{tan2006regression}, \textcite{ogburn2015doubly}, and others; useful overviews and discussions are given by \textcite{hernan2006instruments}, \textcite{imbens2014instrumental} (with discussion), and \textcite{baiocchi2014instrumental}.

\begin{assumption}[Consistency] \label{ass:4consistency}
$A=A^Z$ and $Y=Y^A$ with probability one.
\end{assumption}
\begin{assumption}[Positivity] \label{ass:4positivity}
$(z,\mathbf{x}) \in \text{supp}(Z,\mathbf{X})$ if $\mathbf{x} \in \text{supp}(\mathbf{X})$.
\end{assumption}
\begin{assumption}[Unconfoundedness of $Z$] \label{ass:4unconfound}
$(Y^z,A^z) \perp\!\!\!\perp Z \mid \mathbf{X}$.
\end{assumption}
\begin{assumption}[Exclusion Restriction] \label{ass:4exclusion}
$Y^{za}=Y^a$ with probability one.
\end{assumption}

Consistency means potential treatments $A^z$ and outcomes $Y^a$ are uniquely defined by a unit's own instrument and treatment levels, respectively, and not by others' levels (i.e., no interference), and also not by the way the instrument or treatment are administered (i.e., no different versions). Positivity says that the instrument is not deterministic, in the sense that every unit has some chance of receiving each level of the instrument, regardless of covariates. Unconfoundedness says that the instrument is essentially randomized once we condition on covariates, i.e., that it is unrelated to potential outcomes $Y^z=Y^{zA^z}$ and treatments $A^z$ under different instrument values $Z=z$. The exclusion restriction says that the instrument only affects outcomes through treatment. Assumptions \ref{ass:4consistency}--\ref{ass:4exclusion} can hold by design in trials where the instrument is externally randomized by investigators, but in observational studies these assumptions are typically untestable and require justification based on subject matter.

Finally we also employ the following regularity conditions on the latent threshold distribution and LIV curve.

\begin{assumption}[Instrumentation] \label{ass:4instrument}
$\inf_{t} p(t \mid \mathbf{v}) > 0$.
\end{assumption}
\begin{assumption}[Continuity] \label{ass:4continuity}
$T$ is continuously distributed and $\gamma(t,\mathbf{v})$ is continuous in $t$.
\end{assumption}

Instrumentation means there are at least some units whose treatment would be affected by the value of the instrument (i.e., some who would take treatment when the instrument passes $Z=t$); for now we leave the set over which the infimum is taken ambiguous. 
{
\begin{remark}
A caveat is in order regarding Assumptions 3 and 6: positivity and instrumentation are particularly strong requirements for continuous instruments. Positivity may be violated, for example, if some units simply have no chance at receiving IV levels far away from what they actually received; e.g., in the analysis in Section 5, it may be implausible that some patients who actually live very close to high-level ICUs could ever live very far from them. Similarly, instrumentation can be violated with weak instruments, since moving the IV past certain levels simply may not affect any units' treatment status. Positivity violations can be partially ameliorated by an appropriate choice of weight function (as discussed in Section 3.1), and instrumentation violations by an appropriate choice of the set $\mathcal{T}$ on which the LIV curve is to be estimated; however neither is a panacea. See \textcite{petersen2010diagnosing, westreich2010invited, stock2012survey, baiocchi2014instrumental} and others for relevant discussion of these issues in related settings. Extensions dealing with violations of these assumptions are not considered here. 
\end{remark} }

The following theorem indicates that the LIV curve can be identified with observed data, under the above assumptions.

\begin{theorem}
\label{thm:liv_id}
Suppose Assumption \ref{ass:4thresh} holds. Let $\mathcal{T} \subset \text{supp}(Z)$ denote a compact set on which we wish to identify $\gamma(t,\mathbf{v})$. If Assumptions \ref{ass:4consistency}--\ref{ass:4exclusion} hold for all $z \in \mathcal{T}$ and Assumptions \ref{ass:4instrument}--\ref{ass:4continuity} hold for all $t \in \mathcal{T}$, then the LIV curve is identified for any $t \in \mathcal{T}$ by
\begin{equation}
\label{eq:liv_id}

\gamma(t,\mathbf{v}) =\frac{ \frac{\partial}{\partial z} \mathbb{E}\{ \mathbb{E}(Y \mid \mathbf{X},Z=z) \mid \mathbf{V}=\mathbf{v}\}}{ \frac{\partial}{\partial z} \mathbb{E}\{ \mathbb{E}(A \mid \mathbf{X},Z=z) \mid \mathbf{V}=\mathbf{v}\} } \bigg|_{z=t}.
\end{equation}
\end{theorem}

A proof of Theorem \ref{thm:liv_id} is given in the Supplementary Materials; the logic follows as in more standard settings where $Z$ is discrete and $\mathbf{V}=\mathbf{X}$. Importantly, the LIV curve can only be identified on subsets of $\text{supp}(Z)$; thus as in the binary instrument setting, we cannot identify effects for never-takers or always-takers with $T=\pm\infty$. {As discussed in more detail in the next section, the ratio-of-derivatives structure of the LIV curve makes its study particularly interesting from a theoretical perspective, especially relative to the LATE parameter studied by \textcite{ogburn2015doubly} and others.}

\begin{remark}
From this point forward,  $\gamma(t,\mathbf{v})$ will denote the observed data expression in \eqref{eq:liv_id}, which represents the causal effect given in \eqref{eq:liv} under Assumptions \ref{ass:4consistency}--\ref{ass:4continuity} as described in Theorem \ref{thm:liv_id}. Of course, if the conditions of Theorem \ref{thm:liv_id} do not hold then \eqref{eq:liv_id} may represent something other than the aforementioned causal effect. For example, if only Assumptions 2--4 hold, then we can only think of the instrument as an unconfounded continuous exposure (or dose), and $\gamma(t,\mathbf{v})$ would represent the ratio of derivatives of the dose-response curves $\mathbb{E}(Y^z \mid \mathbf{V}=\mathbf{v})$ and $\mathbb{E}(A^z \mid \mathbf{V}=\mathbf{v})$.
\end{remark}

\section{Main Results}

In this section we develop semiparametric theory for models of the LIV curve defined in \eqref{eq:liv_id}, use this theory to develop novel estimators (including inverse-probability-weighted, regression, and doubly robust estimators), describe asymptotic properties, and finally present cross-validation methods for model selection in high-dimensional settings.

\subsection{Semiparametric Theory}

Suppose we have a parametric model for the LIV curve, which we write as $\gamma(t,\mathbf{v};\boldsymbol\psi)$ for some finite-dimensional $\boldsymbol\psi \in \mathbb{R}^q$. Importantly, we do not assume this model is necessarily correct, and instead follow \textcite{neugebauer2007nonparametric}, \textcite{rosenblum2010targeted}, and others in using a working model approach, by formulating our estimand as the projection of the true curve $\gamma(t,\mathbf{v})$ onto the posed working model. Specifically, we use the weighted least squares projection given by
\begin{equation}
\boldsymbol\psi_0 = \argmin_{\boldsymbol\psi \in \mathbb{R}^q} \ \mathbb{E}\Big[ w(T,\mathbf{V}) \{ \gamma(T,\mathbf{V}) - \gamma(T,\mathbf{V};\boldsymbol\psi) \}^2 \Big] ,
\label{eq:projection}
\end{equation}
where $w(t,\mathbf{v})$ is some user-specified weight function. {We use the $L_2$ loss-based projection for its convenience and familiarity, but it would be worthwhile to develop results for other loss functions in future work.} 

The projection approach warrants some discussion. Whether to use a model like $\gamma(t,\mathbf{v};\boldsymbol\psi)$ only for projections, or to assume it is actually correct, can be viewed as a bias-variance trade-off. If the model happens to be correct, then both the projection and model-based approaches will yield valid estimates of the true function, though projection estimators will not generally be fully efficient (depending on the choice of weight function). However, if the posited model is incorrect, then the model-based approach is technically no longer valid and can be difficult to interpret, since it may not correspond to an interpretable projection (e.g., least-squares); on the other hand, the projection approach is still well-defined and represents a best-fitting summary measure. 

The projection approach essentially formalizes how models are often viewed as approximations in practice. However, we note that there are some distinctions between using approximations in causal versus predictive settings. In causal settings one often cares about understanding mechanisms or making beneficial treatment recommendations; when the model is incorrect, different projections can vary in their usefulness in these respects. And the costs of an inaccurate projection may be quite different for these causal goals (either smaller or larger), relative to usual predictive goals of minimizing prediction error. 

{Another important issue in using projections is choosing the weight $w(t,\mathbf{v})$. Ideally there will be some subject matter justification for learning about specific parts of the effect curve. If not, in theory one could use a uniform weight that assigns mass equally across the support, but in practice this could lead to poor efficiency. We have found that weights based on the instrument density, e.g., $w(t,\mathbf{v})=p(Z=t, \mathbf{v})$ or $w(t,\mathbf{v})=w(t)=p(Z=t)$, work well when no particular weight function is preferred based on substantive concerns. Of course if $\gamma$ is correctly specified, all weights yield consistent estimators that only vary in efficiency.}

Note that the projection parameter in \eqref{eq:projection} depends on the distribution of the latent threshold $T$. Although this threshold is not observed directly, its distribution is identified in the observed data (under Assumptions \ref{ass:4mono}--\ref{ass:4continuity}). For example when $t \in \text{supp}(Z)$ we have
\begin{equation}

p(t \mid \mathbf{v}) =  \frac{\partial}{\partial z} \mathbb{E}\{\mathbb{E}(A \mid \mathbf{X}, Z=z) \mid \mathbf{V}=\mathbf{v}\} \Bigm|_{z=t}
\end{equation}
Importantly, the above expression for the threshold density equals the denominator of  $\gamma(t,\mathbf{v})$ given in Theorem \ref{thm:liv_id}. 

After characterizing the parameter of interest in terms of observed data as in \eqref{eq:projection}, based on the expression in \eqref{eq:liv_id}, it is possible to estimate it using any number of approaches, such as parametric or nonparametric maximum likelihood, or Bayesian methods. In our setting, however, semiparametric approaches have a number of important advantages. First, they can incorporate information about the instrument mechanism, which may be better understood or easier to model than the outcome and treatment mechanisms (which is what a likelihood-based approach would rely on modeling). Second, they allow for double robustness, which means consistent estimation of $\boldsymbol\psi$ is possible as long as either the instrument mechanism or the treatment/outcome mechanisms are correctly modeled (not necessarily all three, so either the instrument or the treatment and outcome models can be misspecified). And third, semiparametric doubly robust approaches allow for fast root-n rates of convergence for the parameter of interest $\boldsymbol\psi$, even when nuisance functions are estimated at slower rates, e.g., using flexible data-adaptive or machine learning methods. Thus these estimators are less sensitive to the curse of dimensionality; this phenomenon was noted recently for example by \textcite{van2014higher} and \textcite{chernozhukov2016double}, who refer to it as orthogonality. Nonetheless our results also lead to novel inverse-weighting and regression-based estimators; the latter in particular might be preferable in small samples.

A crucial aspect of developing semiparametric theory and corresponding estimators for a given problem involves characterizing the possible influence functions, and in particular finding the efficient influence function. Many details on semiparametric theory are available elsewhere \autocite{bickel1993efficient, van2003unified, tsiatis2006semiparametric, kennedy2016semiparametric}, so we give only a brief review here. Any regular asymptotically linear estimator minus its target parameter can be expressed as the empirical average of its so-called influence function plus an $o_\mathbb{P}(1/\sqrt{n})$ error term. Viewed as elements of a Hilbert space of mean-zero finite-variance functions equipped with covariance norm, the influence functions under a given model lie in the orthogonal complement of the nuisance tangent space. The efficient influence function can then be defined as the influence function with smallest variance, the projection of any influence function onto the tangent space of scores, or as a particular pathwise derivative. The efficient influence function is especially important in practice because its variance is the semiparametric efficiency bound (thus providing a benchmark for efficient estimation), and because it can be used to construct estimators that are doubly robust and potentially semiparametric efficient.

A major challenge in deriving semiparametric theory for the projection parameter in \eqref{eq:projection} is its complexity; namely, it is a weighted projection of a ratio of derivatives of regression functions that are partially marginalized. {We conjecture that this may be why a doubly robust estimator for the LIV curve has yet to appear in the literature. In general, such complex structure would yield a complicated efficient influence function involving derivatives of regression functions, making corresponding estimators very difficult to compute. A major contribution of our work is an expression for the efficient influence function that only involves derivatives based on the known (and generally more simple) model and weight functions, rather than unknown complex and high-dimensional nuisance functions. Although still complex, our formulation yields efficient estimators that are easier to construct in practice.} The next theorem gives the efficient influence function for the parameters of the LIV curve projection.

\begin{theorem} \label{thm:4eif}
Suppose the weight function $w(t,\mathbf{v})$ is continuously differentiable in $t$ and satisfies $w(t,\mathbf{v})=0$ for $t \notin \text{int}(\mathcal{T})$, with the set $\mathcal{T} \subset \text{supp}(Z)$ defined as in Theorem \ref{thm:liv_id}. Also assume that partial derivatives (with respect to $\boldsymbol\psi$ and $t$) of the working model $\gamma(t,\mathbf{v};\boldsymbol\psi)$ exist and are continuous. Then, under a nonparametric model, the efficient influence function for $\boldsymbol\psi$ defined in  \eqref{eq:projection} is proportional to
\begin{align} \label{eq:eif}
\mathbf{o}ldsymbol\mathbb{V}phi(\mathbf{O};\boldsymbol\psi,\mathbf{o}ldsymbol\eta) &=  \int_\mathcal{T} \Big\{ \mathbf{g}_1(t,\mathbf{V};\boldsymbol\psi) \mathbb{E}(A \mid \mathbf{X},Z=t) - \mathbf{g}_2(t,\mathbf{V};\boldsymbol\psi) \mathbb{E}(Y \mid \mathbf{X},Z=t) \Big\} \ dt  \\
& \hspace{.4in} + \mathbf{g}_1(Z,\mathbf{V}; \boldsymbol\psi) \left\{ \frac{A-\mathbb{E}(A \mid \mathbf{X},Z)}{p(Z \mid \mathbf{X})} \right\} - \mathbf{g}_2(Z,\mathbf{V};\boldsymbol\psi) \left\{ \frac{Y - \mathbb{E}(Y \mid \mathbf{X},Z)}{p(Z \mid \mathbf{X})} \right\} \nonumber  
\end{align}
where $\mathbf{o}ldsymbol\eta=(\pi,\lambda,\mu)$ denotes the nuisance functions defined in Section \ref{4data}, and $\mathbf{g}_1$ and $\mathbf{g}_2$ are the $(q \times 1)$ vectors
\begin{align*}
\mathbf{g}_1(z,\mathbf{v};\boldsymbol\psi) &= \frac{\partial}{\partial t} \left\{ \frac{\partial}{\partial \boldsymbol\psi^*} \gamma(t,\mathbf{v}; \boldsymbol\psi^*) \Big|_{\boldsymbol\psi^*=\boldsymbol\psi} w(t,\mathbf{v}) \gamma(t,\mathbf{v};\boldsymbol\psi) \right\}\Big|_{t=z} \\
\mathbf{g}_2(z,\mathbf{v};\boldsymbol\psi) &= \frac{\partial}{\partial t} \left\{ \frac{\partial}{\partial \boldsymbol\psi^*} \gamma(t,\mathbf{v}; \boldsymbol\psi^*)  \Big|_{\boldsymbol\psi^*=\boldsymbol\psi} w(t,\mathbf{v})  \right\} \Big|_{t=z} .
\end{align*}
\end{theorem}

A proof of Theorem \ref{thm:4eif} is given in the Supplementary Materials. The first step is to define $\boldsymbol\psi$ as the zero of a moment condition coming from the derivative of \eqref{eq:projection}. Then, to deal with the crucial difficulty that this moment condition involves complex derivatives of partially marginalized treatment/outcome regression functions as in \eqref{eq:liv_id}, we use integration by parts to instead transfer the derivatives  to the known model $\gamma(t,\mathbf{v};\boldsymbol\psi)$ and weight function $w(t,\mathbf{v})$. This allows us to write the moment condition as a somewhat more standard causal quantity, essentially equaling sums of the mean treatment/outcome under a stochastic intervention on the instrument, where the conditional instrument density is replaced with a uniform distribution on $\mathcal{T}$.  At this point one can proceed by showing that the function $\mathbf{o}ldsymbol\mathbb{V}phi$ is the canonical gradient of the pathwise derivative of $\boldsymbol\psi$. 

This explains why the influence function in Theorem \ref{thm:4eif} takes the form it does: it resembles that of a (uniform) stochastic intervention effect, as in for example \textcite{diaz2012population}, with the non-standard form of the functions $\mathbf{g}_1$ and $\mathbf{g}_2$ coming from the integration by parts trick.  The influence function as a whole can be viewed as the ``derivative term'' in a generalized von Mises expansion of the moment condition, as in for example \textcite{van2014higher}; alternatively, as in other causal inference and missing data problems, it can be viewed as consisting of an inverse-probability-weighted term (the added terms in the second line of \eqref{eq:eif}) plus an augmentation term (the first line and subtracted terms in the second line of \eqref{eq:eif}).

Importantly, using integration by parts to transfer derivatives from $P$-dependent quantities also means the influence function can be evaluated without analytical differentiation of the regression functions, which makes it much more practical for constructing and implementing estimators, as discussed in more detail in the next subsection. The user-specified weight function is required to vanish outside the interior of the set $\mathcal{T} \subset \text{supp}(Z)$ since the LIV curve $\gamma(t,\mathbf{v})$ is not identified outside of $\mathcal{T}$ as discussed in Theorem \ref{thm:liv_id}. 

\subsection{Proposed Method}

Once we have derived the efficient influence function, we can use it to construct estimators that have numerous advantageous properties. A standard approach is to solve an estimating equation based on an estimated version of the efficient influence function; specifically we can use $\mathbf{o}ldsymbol\mathbb{V}phi$ as an estimating function, with unknown nuisance functions replaced with estimates. 

Thus our proposed estimator for a given working model $\gamma(t,\mathbf{v};\boldsymbol\psi)$ is given by $\mathbf{o}ldsymbol{\hat\psi}$, defined as the solution in $\boldsymbol\psi$ to the estimating equation
\begin{equation} \label{eq:4est}
 \mathbb{P}_n\{ \mathbf{o}ldsymbol\mathbb{V}phi(\mathbf{O}; \boldsymbol\psi, \mathbf{o}ldsymbol{\hat\eta}) \} = \mathbf{0} , 
\end{equation}
where $\mathbf{o}ldsymbol{\hat\eta}=(\hat\pi,\hat\lambda,\hat\mu)$ are estimated versions of the three nuisance functions $\mathbf{o}ldsymbol{\eta}=(\pi,\lambda,\mu)$. An inverse-probability-weighted estimator can be obtained by using $\hat\lambda=\hat\mu=0$, and a regression estimator can be obtained by using $\hat\pi=\infty$. Another option for constructing estimators based on influence functions is targeted minimum loss-based methodology \autocite{van2006targeted}, which yields plug-in estimators that respect the bounds of the parameter space. In our setting, our proposed estimating equation approach will also respect any such bounds, as long as the chosen working model does; it is also relatively straightforward to implement. 
{
\begin{remark}
For notational simplicity, the estimator proposed above uses nuisance estimates constructed from the entire sample (with asymptotic results in Section 3.3 relying on empirical process conditions). In Section 4 of the Supplementary Materials we present sample-splitting estimation, in the same spirit as \textcite{robins2008higher}, \textcite{zheng2010asymptotic},  \textcite{robins2013new}, \textcite{chernozhukov2016double}, and others, which does not require empirical process conditions and allows arbitrarily complex and adaptive nuisance estimators, e.g., random forests.
\end{remark} }

First consider the simple case where $\mathbf{V}=\emptyset$ (i.e., effect modification is not of interest), and the LIV curve $\gamma(t)$ is projected onto a constant $\gamma(t;\boldsymbol\psi)=\psi$. Here the target estimand $\psi$ is a simple weighted average of $\gamma(t)$, of the form
$$ \psi = \int_\mathcal{T} w^*(t) \gamma(t) \ dt $$
with weight $w^*(t)= w(t) p(t) / \int_\mathcal{T} w(t) p(t) \ dt$, {where $w(t)=w(t,\mathbf{v})$ since $\mathbf{V}$ is empty}. Note that in general the parameter $\boldsymbol\psi$ indexes a model for the threshold-dependent LIV curve $\gamma(t,\mathbf{v})$, but here it is a weighted average over $T$ because we have used a simple projection and are agnostic about whether the constant effect model is correct. In this case the quantity $\psi$ can also be viewed as the mean treatment effect (among compliers) in a population where the density of the latent threshold $T$ among compliers equals $w^*(t)$. Solving \eqref{eq:4est} leads to the ratio estimator
\begin{equation} \label{eq:wavg}
 \hat\psi = \frac{ \int_\mathcal{T} w'(t) \hat{m}(t) \ dt + \mathbb{P}_n \left\{ w'(Z) \frac{Y-\hat\mu(\mathbf{X},Z)}{\hat\pi(Z \mid \mathbf{X})} \right\}  } { \int_\mathcal{T} w'(t) \hat\ell(t) \ dt + \mathbb{P}_n \left\{ w'(Z)  \frac{A-\hat\lambda(\mathbf{X},Z)}{\hat\pi(Z \mid \mathbf{X})} \right\} } 
\end{equation}
where $\hat{m}(t)=\mathbb{P}_n\{\hat\mu(\mathbf{X},t) \}$ and $\hat\ell(t)=\mathbb{P}_n \{ \hat\lambda(\mathbf{X},t) \}$ are estimates of the marginalized regression functions from Section \ref{4data}. Thus $\hat\psi$ is an adjusted version of the regression-based plug-in estimator $\int_\mathcal{T} w'(t) \hat{m}(t) \ dt /  \int_\mathcal{T} w'(t) \hat\ell(t) \ dt$, where adding inverse-probability-weighted terms to the numerator and denominator is the adjustment required to obtain double robustness. 

More standard IV estimators are often computed with a two-stage least squares approach, where in the first stage the treatment variable is regressed on the instrument (and covariates) and then in the second stage the outcome is regressed on the predicted values from the first stage (and covariates). In fact, the weighted average estimator in \eqref{eq:wavg} can also be constructed with a modified version of such a two-stage least squares approach; this may make it more amenable to practical use. Specifically, the following modified two-stage least squares procedure can be used to compute the weighted average estimator  (using pseudo- instrument, treatment, and outcome $w'(Z)$, $A^*$, and $Y^*$ respectively):

\begin{enumerate}
\item Regress $A^* = \frac{A - \hat\lambda(\mathbf{X},Z)}{\pi(Z \mid \mathbf{X})} + \frac{ \mathbbold{1}\{w'(Z) \neq 0\}}{w'(Z)} \int_\mathcal{T} w'(t) \hat\lambda(\mathbf{X},t) \ dt$ on $w'(Z)$ without an intercept, and obtain predicted values $\hat{A}^*$. 
\item Regress $Y^* = \frac{Y - \hat\mu(\mathbf{X},Z)}{\pi(Z \mid \mathbf{X})} + \frac{ \mathbbold{1}\{w'(Z) \neq 0\}}{w'(Z)} \int_\mathcal{T} w'(t) \hat\mu(\mathbf{X},t) \ dt$ on $\hat{A}^*$, without an intercept.
\end{enumerate}
Then the coefficient in front of $\hat{A}^*$ in the second stage equals $\hat\psi$ from \eqref{eq:wavg}.

Closed-form estimators are also available even when effect modification is of interest, as long as we project onto linear models of the form $\gamma(t,\mathbf{v};\boldsymbol\psi)=\mathbf{h}(t,\mathbf{v})^{ \mathrm{\scriptscriptstyle T} } \boldsymbol\psi$, for some known mapping $\mathbf{h}: \mathcal{T} \times \text{supp}(\mathbf{V}) \rightarrow \mathbb{R}^q$. Specifically, in such cases the estimator $\mathbf{o}ldsymbol{\hat\psi}$ defined as the solution to \eqref{eq:4est} is given by
\begin{align*}
\mathbf{o}ldsymbol{\hat\psi} &= \mathbb{P}_n \left[ \mathbf{g}_1^*(Z,\mathbf{V}) \left\{ \frac{A-\hat\lambda(\mathbf{X},Z)}{\hat\pi(Z \mid \mathbf{X})} \right\} + \int_\mathcal{T} \mathbf{g}_1^*(t,\mathbf{V}) \hat\lambda(\mathbf{X},t) \ dt \right]^{-1} \\
& \hspace{.4in} \times \mathbb{P}_n \left[ \mathbf{g}_2(Z,\mathbf{V}) \left\{ \frac{Y-\hat\mu(\mathbf{X},Z)}{\hat\pi(Z \mid \mathbf{X})} \right\} + \int_\mathcal{T} \mathbf{g}_2(t,\mathbf{V}) \hat\mu(\mathbf{X},t) \ dt \right]
\end{align*}
where $\mathbf{g}_1^*(z,\mathbf{v}) = \frac{\partial}{\partial t} \{ \mathbf{h}(t,\mathbf{v}) w(t,\mathbf{v}) \mathbf{h}(t,\mathbf{v})^{ \mathrm{\scriptscriptstyle T} }\}|_{t=z}$, and $\mathbf{g}_2(z,\mathbf{v})=\frac{\partial}{\partial t} \{ \mathbf{h}(t,\mathbf{v}) w(t,\mathbf{v})\}|_{t=z}$ is as defined in Theorem \ref{thm:4eif}. Closed-form expressions will typically not be available for estimators in general non-linear models; however, since such estimators are still defined as estimating equation-based Z-estimators, they can be computed with standard software (for example, one could use the \verb+optim+ function in R). Variance estimation and confidence interval construction will be discussed in the next section.

\subsection{Asymptotic Theory}

In this section we discuss the asymptotic properties of our proposed estimation approach. In particular we show that our estimator is doubly robust, and that if the nuisance functions are estimated well enough it is asymptotically normal and efficient. Further, asymptotic normality and efficiency are possible even after flexible machine learning-based covariate adjustment. (Our results equally apply to estimators that only solve the efficient influence function estimating equation asymptotically, up to order $o_\mathbb{P}(1/\sqrt{n})$, such as targeted minimum loss-based estimators.) {An analogous result holds for a sample-splitting version of the proposed estimator, as detailed in Section 6 of the Supplementary Materials.}

\begin{theorem} \label{thm:4asymptotics}
Assume that:
\begin{enumerate}
\item $(\mathbf{o}ldsymbol{\hat\psi},\mathbf{o}ldsymbol{\hat\eta}) \stackrel{p}{\rightarrow} (\boldsymbol\psi_0,\mathbf{o}ldsymbol{\overline\eta})$, where $\mathbf{o}ldsymbol{\overline\eta}=(\overline\pi,\overline\lambda,\overline\mu)$ with either $\overline\pi=\pi_0$ or $(\overline\lambda,\overline\mu)=(\lambda_0,\mu_0)$.

\item The sequence of functions $\mathbf{o}ldsymbol{\hat\mathbb{V}phi}_n = \mathbf{o}ldsymbol\mathbb{V}phi(\cdot;\mathbf{o}ldsymbol{\hat\psi},\mathbf{o}ldsymbol{\hat\eta})$ and its limit $\mathbf{o}ldsymbol\mathbb{V}phi_0=\mathbf{o}ldsymbol\mathbb{V}phi(\cdot;\mathbf{o}ldsymbol{\psi_0},\mathbf{o}ldsymbol{\overline\eta})$ are contained in a Donsker class with $||\mathbf{o}ldsymbol{\hat\mathbb{V}phi}_n - \mathbf{o}ldsymbol\mathbb{V}phi_0||_2=o_\mathbb{P}(1)$.
\item The map $\boldsymbol\psi \rightarrow \mathbb{E}\{ \mathbf{o}ldsymbol\mathbb{V}phi(\mathbf{O};\boldsymbol\psi,\mathbf{o}ldsymbol\eta) \}$ is differentiable at $\boldsymbol\psi_0$ uniformly in $\mathbf{o}ldsymbol\eta$ (around $\mathbf{o}ldsymbol{\overline\eta}$), with invertible derivative matrix $\mathbf{D}(\boldsymbol\psi_0,\mathbf{o}ldsymbol\eta) \rightarrow \mathbf{D}(\boldsymbol\psi_0,\mathbf{o}ldsymbol{\overline\eta}) \equiv \mathbf{D}_0$.
\end{enumerate}
Then the proposed estimator is consistent with rate of convergence
$$ || \mathbf{o}ldsymbol{\hat\psi} - \boldsymbol\psi_0 || = O_\mathbb{P}\left\{ 1/\sqrt{n} + ||\hat\pi - \pi_0||_2 \Big( ||\hat\lambda - \lambda_0 ||_2 + ||\hat\mu - \mu_0 ||_2 \Big) \right\} . $$
Suppose further that:
\begin{enumerate}
\item[(d)] $||\hat\pi - \pi_0||_2 ( ||\hat\lambda - \lambda_0 ||_2 + ||\hat\mu - \mu_0 ||_2 )=o_\mathbb{P}(1/\sqrt{n})$.
\end{enumerate}
Then the proposed estimator is asymptotically normal with
$$ \sqrt{n}(\mathbf{o}ldsymbol{\hat\psi} - \boldsymbol\psi_0) \rightsquigarrow N\Big(\mathbf{0}, \mathbb{E}[\{ \mathbf{D}_0^{-1} \mathbf{o}ldsymbol\mathbb{V}phi(\mathbf{O};\boldsymbol\psi_0,\mathbf{o}ldsymbol\eta_0) \}^{\otimes 2}] \Big) , $$
and thus semiparametric efficient.
\end{theorem}

A proof of Theorem \ref{thm:4asymptotics} is given in the Supplementary Materials; it follows from standard Z-estimator theory and empirical process results \autocite{van1996weak,  van2002semiparametric}. The first condition indicates the double robustness of our approach, since some of the nuisance estimators $\mathbf{o}ldsymbol{\hat\eta}=(\hat\pi,\hat\lambda,\hat\mu)$ can be misspecified. Specifically, as long as either $\hat\pi$ or $(\hat\lambda,\hat\mu)$ is consistent, then the estimator $\mathbf{o}ldsymbol{\hat\psi}$ will be as well. This gives analysts two chances at consistency, and is particularly important in the IV setting since it can be easier to model the instrument density $\pi$ than the two regression functions $(\lambda,\mu)$. 

Conditions (b)--(c) of Theorem \ref{thm:4asymptotics} are standard regularity conditions for M- and Z-estimators  \autocite{van1996weak, van2000asymptotic, van2002semiparametric}. {Condition (b) restricts the flexibility of the nuisance estimators (and their limits), but Donsker classes still cover many complex functions.} For example, parametric Lipschitz functions are Donsker, but so are many more complicated function types such as infinite-dimensional smooth functions with bounded partial derivatives, VC classes, Sobolev classes, and functions with bounded uniform sectional variation, as well as convex combinations and Lipschitz transformations of any these classes. More discussion and examples can be found in Sections 2.6--2.7 of \textcite{van1996weak} and Examples 19.6--19.12 of \textcite{van2000asymptotic}, as well as in \textcite{kennedy2016semiparametric}. Condition (b) is important because it means we do not have to rely on restrictive parametric models to estimate the potentially complicated and high-dimensional nuisance functions $(\pi,\lambda,\mu)$, and can instead use more flexible data-adaptive methods. Condition (b) can also be weakened in various ways. For example, the Donsker condition really only needs to hold in a shrinking neighborhood of $(\boldsymbol\psi_0,\mathbf{o}ldsymbol{\overline\eta})$, or with high probability as $n \rightarrow \infty$; alternatively we could formulate Condition (b) in terms of weaker entropy or bracketing conditions. {Importantly, in Section 6 of the Supplementary Materials, we discuss how the sample-splitting estimator mentioned earlier (and defined in Section 4 of the Supplementary Materials) can do away with empirical process conditions entirely, allowing arbitrarily flexible nuisance estimators such as random forests}. The differentiability in Condition (c) is standard and required to use a delta method-type result (note that the influence function need not be differentiable itself, only its expectation).

Under Conditions (a)--(c) of Theorem \ref{thm:4asymptotics}, the proposed estimator is consistent with rate of convergence given by $1/\sqrt{n} + ||\hat\pi - \pi_0||_2 ( ||\hat\lambda - \lambda_0 ||_2 + ||\hat\mu - \mu_0 ||_2 )$. Again the double robustness is apparent since consistency (i.e., $|| \mathbf{o}ldsymbol{\hat\psi} - \boldsymbol\psi_0 || = o_\mathbb{P}(1)$) follows as long as either the instrument density is consistently estimated, i.e., $||\hat\pi - \pi_0||_2 = o_\mathbb{P}(1)$, or the treatment/outcome regressions are, i.e., $(||\hat\lambda - \lambda_0 ||_2 + ||\hat\mu - \mu_0 ||_2)=o_\mathbb{P}(1)$. 
Importantly, however, the result also shows how double robustness is useful even apart from giving two chances at consistency; in particular, if we estimate the regression functions $(\lambda,\mu)$ at slower rates, double robustness gives us a chance to obtain faster rates for $\mathbf{o}ldsymbol{\hat\psi}$ by consistently estimating $\pi$, and vice versa. 

For example, if Condition (d) holds so that $||\hat\pi - \pi_0||_2 ( ||\hat\lambda - \lambda_0 ||_2 + ||\hat\mu - \mu_0 ||_2 )=o_\mathbb{P}(1/\sqrt{n})$, i.e., effects of nuisance estimation are asymptotically negligible, then the estimator $\mathbf{o}ldsymbol{\hat\psi}$ is root-n consistent, asymptotically normal, and semiparametric efficient. Note Condition (d) can hold even if the nuisance functions are estimated at slower than parametric root-n rates, so that efficient estimation and valid inference is possible for $\boldsymbol\psi$ even if we use machine learning-based covariate adjustment, via flexible estimation of the nuisance functions $(\pi,\lambda,\mu)$. For example, if the nuisance functions $(\pi,\lambda,\mu)$ are all estimated at faster than $n^{1/4}$ rates, so that $||\hat\pi - \pi_0||_2 = ||\hat\lambda - \lambda_0 ||_2 = ||\hat\mu - \mu_0 ||_2 =o_\mathbb{P}(n^{-1/4})$, then Condition (d) holds since $o_\mathbb{P}(n^{-1/4})o_\mathbb{P}(n^{-1/4})=o_\mathbb{P}(1/\sqrt{n})$. {Thus in this case the rate for $\hat{\boldsymbol\psi}$ is faster than that of the nuisance estimators  (here, it is the square of the slower nuisance rates).} Such $n^{1/4}$ rates are possible in various flexible models; for instance, under some conditions \autocite{horowitz2009semiparametric} generalized additive model estimators can obtain rates of the form $O_\mathbb{P}(n^{-2/5})$, which is $o_\mathbb{P}(n^{-1/4})$ since 
$$R_n=O_\mathbb{P}(n^{-2/5}) \implies n^{1/4}R_n = n^{-3/20} n^{2/5} R_n=O_\mathbb{P}(n^{-3/20})=o_\mathbb{P}(1). $$ 
Condition (d) can also hold if one of $\pi$ or $(\lambda,\mu)$ is estimated with a correctly specified parametric model and the other is merely estimated consistently. 

If Condition (d) holds, confidence intervals can be constructed with the bootstrap, or using a direct estimate of the asymptotic variance given in Theorem \ref{thm:4asymptotics}, such as

$$  \mathbb{P}_n[ \{ \mathbf{\hat{D}}^{-1}\mathbf{o}ldsymbol\mathbb{V}phi(\mathbf{O};\mathbf{o}ldsymbol{\hat\psi},\mathbf{o}ldsymbol{\hat\eta})\}^{\otimes 2}] $$
where $\mathbf{\hat{D}}= \mathbb{P}_n\{ \partial \mathbf{o}ldsymbol\mathbb{V}phi(\mathbf{O};\mathbf{o}ldsymbol{\psi},\mathbf{o}ldsymbol{\hat\eta}) / \partial \boldsymbol\psi^{ \mathrm{\scriptscriptstyle T} }\}|_{\boldsymbol\psi=\mathbf{o}ldsymbol{\hat\psi}}$ is an estimate of the derivative matrix from Condition (c) of Theorem \ref{thm:4asymptotics}. {For completeness}, we note that if parametric models are used to estimate all three nuisance functions (and either $\pi$ or $(\lambda,\mu)$ are correctly modeled, but not both) then the bootstrap would still be valid even though Condition (d) fails, since then the contribution from nuisance estimation is asymptotically linear (an analytic expression could also be derived, since $\mathbf{o}ldsymbol{\hat\psi}$ and the estimated nuisance parameters solve a large system of estimating equations). {However, as we note in the Introduction, parametric nuisance models are likely to be misspecified, except for $\pi$ if $Z$ is externally randomized.} Inference is somewhat more delicate in a truly doubly robust but nonparametric setting, where one nuisance estimator can be misspecified but $(\hat\pi,\hat\lambda,\hat\mu)$ are estimated flexibly; we leave this to future work. 

\subsection{Model Selection}

To this point we have presumed that we have an a priori model $\gamma(t,\mathbf{v};\boldsymbol\psi)$, which either represents the truth or a low-dimensional projection. When such a priori models are not available, we might instead aim to learn the form of $\gamma(t,\mathbf{v})$ from data. Thus in this section we propose a doubly robust cross-validation approach for model selection. Model selection is an important issue in causal inference in general, but this is especially the case for the LIV curve, since the latent threshold $T$ is continuous; thus saturated parametric models are not possible, even when effect modification is not of interest (i.e., $\mathbf{V}=\emptyset$). In this section we derive the efficient influence function for the risk of a given candidate estimator, and show how it can be used as a doubly robust loss function in the cross-validation framework developed by \textcite{van2003cross}. {In particular, our cross-validation approach could be used to select among working models whose complexity increases with sample size, yielding doubly robust yet nonparametric estimators of the local IV curve, in the same spirit as \textcite{robins2001inference} and \textcite{kennedy2017nonparametric}.}

If we knew the true LIV curve and true distribution $P$, we could evaluate the performance of a given estimator $\hat\gamma_k$ by computing the mean squared error risk
$R^*(\hat\gamma_k) = \int_\mathcal{V} \int_\mathcal{T} w(t,\mathbf{v}) \{ \gamma(t,\mathbf{v}) - \hat\gamma_k(t,\mathbf{v}) \}^2 \ dP(t,\mathbf{v})$. 
Alternatively, if we only wanted to compare or rank a set of candidate estimators $\{ \hat\gamma_k : k \in \mathcal{K}\}$, we could use the pseudo-risk
\begin{equation} \label{eq:4risk}
R(\hat\gamma_k) = \int_\mathcal{V} \int_\mathcal{T} w(t,\mathbf{v}) \Big\{ \hat\gamma_k(t,\mathbf{v})^2 - 2\gamma(t,\mathbf{v}) \hat\gamma_k(t,\mathbf{v})  \Big\} \ dP(t,\mathbf{v}) ,
\end{equation}
since $R(\hat\gamma_k)=R^*(\hat\gamma_k) - \mathbb{E}\{w(T,\mathbf{V}) \gamma(T,\mathbf{V})^2\}$ is simply a shifted version of the mean squared error $R^*(\hat\gamma_k)$, and the shift does not depend on the candidate estimator $\hat\gamma_k$. In standard cross-validation it is possible to estimate risk unbiasedly, without worrying about nuisance function estimation; in contrast, in our setting the risk parameter $R(\hat\gamma_k)$ depends on complex nuisance functions via the curve $\gamma(t,\mathbf{v})$ and the distribution of the threshold $T$. Thus estimation of the risk $R(\hat\gamma_k)$ itself requires nuisance estimation, and in fact we can treat $R(\hat\gamma_k)$ as a parameter in its own right, for which we can develop semiparametric theory and estimators. Thus in the next theorem we give the efficient influence function for the risk $R(\gamma_k)$  for a given fixed candidate $\gamma_k$, and go on to show how to use this efficient influence function as a doubly robust loss function for cross-validation-based model selection. 

\begin{theorem} \label{thm:4eif2}
Consider the same setting and assumptions as in Theorem \ref{thm:4eif}. Under a nonparametric model, the efficient influence function for the risk $R(\gamma_k)$ defined in \eqref{eq:4risk} for a fixed candidate $\gamma_k$ is given by $L(\mathbf{O};\gamma_k,\mathbf{o}ldsymbol\eta) - R(\gamma_k)$, for
\begin{align} \label{eq:4loss}
L(\mathbf{O};\gamma_k,\mathbf{o}ldsymbol\eta) &= \int_\mathcal{T} \Big\{ f_1(t,\mathbf{V};\gamma_k) \mathbb{E}(Y \mid \mathbf{X},Z=t) - f_2(t,\mathbf{V};\gamma_k) \mathbb{E}(A \mid \mathbf{X}, Z=t) \Big\} \ dt  \\
& \hspace{.4in} + f_1(Z, \mathbf{V}; \gamma_k) \left\{ \frac{Y - \mathbb{E}(Y \mid \mathbf{X},Z)}{p(Z \mid \mathbf{X})} \right\}  - f_2(Z,\mathbf{V};\gamma_k) \left\{ \frac{A-\mathbb{E}(A \mid \mathbf{X},Z)}{p(Z \mid \mathbf{X})} \right\}  \nonumber 
\end{align}
where $\mathbf{o}ldsymbol\eta=(\pi,\lambda,\mu)$ are the nuisance functions from before, and $f_1$ and $f_2$ are defined as
\begin{align*}
f_1(z,\mathbf{v};\boldsymbol\psi) &= 2\frac{\partial}{\partial t} \Big\{w(t,\mathbf{v}) \gamma_k(t,\mathbf{v}) \Big\}\Big|_{t=z} \\ 
f_2(z,\mathbf{v};\boldsymbol\psi) &= \frac{\partial}{\partial t} \Big\{ w(t,\mathbf{v}) \gamma_k(t,\mathbf{v})^2 \Big\} \Big|_{t=z} . 
\end{align*}
\end{theorem}

A proof of Theorem \ref{thm:4eif2} is given in the Supplementary Materials, and follows similar logic as the proof of Theorem \ref{thm:4eif}. We also show that $L(\mathbf{O};\gamma_k,\mathbf{o}ldsymbol\eta)$ is a doubly robust loss function for the risk $R(\gamma_k)$ in the sense that $\mathbb{E}\{L(\mathbf{O};\gamma_k,\mathbf{o}ldsymbol{\overline\eta})\}=R(\gamma_k)$ for nuisance function $\mathbf{o}ldsymbol{\overline\eta}=(\overline\pi,\overline\lambda,\overline\mu)$ as long as either $\overline\pi=\pi_0$ or $(\overline\lambda,\overline\mu)=(\lambda_0,\mu_0)$, and not necessarily both. Thus we can use $L(\mathbf{O};\gamma_k,\mathbf{o}ldsymbol\eta)$ as a doubly robust estimating function, similar to how we used $\mathbf{o}ldsymbol\mathbb{V}phi(\mathbf{O};\boldsymbol\psi,\mathbf{o}ldsymbol\eta)$ in previous sections. However, since we typically do not have an independent sample to generate candidates $\hat\gamma_k$, we need to generate them from the same sample in which we estimate risk. Thus we can use sample-splitting  to prevent over-fitting.

In particular, we propose using the loss function in \eqref{eq:4loss} for doubly robust model selection following the general approach of \textcite{van2003cross}. This requires some new notation. Let $\mathbf{S}=(S_1,...,S_n)$ denote a random variable independent of the sample that splits the data into training ($S_i=0$) and test ($S_i=1$) sets. For example standard $v$-fold cross-validation arises by allowing the split variable $\mathbf{S}$ to take $v$ different values $\{ \mathbf{S}_1,...,\mathbf{S}_v\}$, each with equal probability $1/v$, where $\sum_i S_{iv}=n/v$ for all $v$ and $\sum_v S_{iv}=1$ for all $i$, so that test sets are all of size $n/v$ and each unit is only used in one test set. Further define $\mathbb{P}_\mathbf{s}^0$ and $\mathbb{P}_\mathbf{s}^1$ as the sub-empirical distributions for the training data $\{i: S_i=0\}$ and test data $\{i: S_i=1\}$, respectively, for a given split $\mathbf{S}=\mathbf{s}$. Therefore, for example, $\mathbf{o}ldsymbol{\hat\eta}(\mathbb{P}_\mathbf{s}^0)$ denotes the nuisance function estimates based only on the training set data, and $\hat\gamma_k(\mathbb{P}_\mathbf{s}^0)$ denotes the LIV curve estimate based only on the training set data (which also depends on the nuisance function estimates constructed from the training data).

The cross-validation selection approach of \textcite{van2003cross} is very similar to standard cross-validation, but incorporates extra steps for nuisance function estimation; it proceeds as follows. For a given split $\mathbf{s}$ and candidate estimator $\hat\gamma_k$, we first estimate the nuisance functions with the training data to obtain $\mathbf{o}ldsymbol{\hat\eta}(\mathbb{P}_\mathbf{s}^0)$, and then estimate the LIV curve with the training data to obtain $\hat\gamma_k(\mathbb{P}_\mathbf{s}^0)$. Then the loss function $L$ can be evaluated for any observation $\mathbf{O}_i$ based on these training estimates, and thus we do so on the test data $\mathbb{P}_\mathbf{s}^1$ and compute the average, given by
$$ \hat{R}_\mathbf{s}(\hat\gamma_k) = \int L\Big\{\mathbf{o}; \hat\gamma_k(\mathbb{P}_\mathbf{s}^0), \mathbf{o}ldsymbol{\hat\eta}(\mathbb{P}_\mathbf{s}^0) \Big\} \ d\mathbb{P}_\mathbf{s}^1(\mathbf{o}) , $$
which we call the estimated risk for candidate $k$ at the current split $\mathbf{s}$. We repeat the above process for each split, average the split-specific risk estimates to get an overall risk estimate for candidate $k$, defined as $\hat{R}(\hat\gamma_k) = \mathbb{E}_\mathbf{S}\{ \hat{R}_\mathbf{S}(\hat\gamma_k) \}$, and finally repeat for each candidate $k \in \mathcal{K}$ and pick the one $\hat{k}$ that yields the smallest overall risk estimate $\hat{k}=\argmin_{k \in \mathcal{K}} \hat{R}(\hat\gamma_k)$. Hence the cross-validation selector can be written as
\begin{equation}
\hat{k} = \argmin_{k \in \mathcal{K}} \ \mathbb{E}_\mathbf{S} \int L\Big\{\mathbf{o}; \hat\gamma_k(\mathbb{P}_\mathbf{S}^0), \mathbf{o}ldsymbol{\hat\eta}(\mathbb{P}_\mathbf{S}^0) \Big\} \ d\mathbb{P}_\mathbf{S}^1(\mathbf{o}) .
\end{equation}
 \textcite{van2003cross} gave conditions under which the risk $\hat{R}(\hat\gamma_{\hat{k}})$ of the above cross-validation selector is asymptotically equivalent to that of an oracle selector given by $\tilde{k} = \argmin_{k \in \mathcal{K}} \ \mathbb{E}_\mathbf{S} \int L\{\mathbf{o}; \hat\gamma_k(\mathbb{P}_\mathbf{S}^0), \mathbf{o}ldsymbol{\overline\eta} \} \ dP(\mathbf{o})$, 
along with corresponding finite-sample bounds. {One important condition is that the number of candidates  does not grow faster than polynomially with sample size, i.e., $|\mathcal{K}| \leq n^c$.} We refer to \textcite{van2003cross} for more details. {Note that semiparametric doubly robust post-selection inference is not considered here, but is an important avenue for future work}.

\section{Simulation Study}

To explore finite sample properties of our methods, we simulated from a model with
\begin{equation*}
\begin{gathered}
(Y^0, \mathbf{X} )\sim N(\mathbf{0},\mathbf{I}_5), \\
Z \mid \mathbf{X}, Y^0 \sim \text{TN}\{1.5 \times \text{sign}(\mathbf{o}ldsymbol\alpha^T \mathbf{X}),4,(-2,2)\}, \\ 
T \mid Z, \mathbf{X}, Y^0 \sim  N(\mathbf{o}ldsymbol\beta^T \mathbf{X}+Y^0,1), \\
A = \mathbbold{1}(Z \geq T), \
Y = Y^0 + A (\psi T),
\end{gathered}
\end{equation*}
where TN$\{\mu,\sigma^2,(l,u)\}$ denotes a truncated normal distribution with support $[l,u]$,  $\mathbf{o}ldsymbol\alpha=(1,1,-1,-1)^{ \mathrm{\scriptscriptstyle T} }$, $\mathbf{o}ldsymbol\beta=(1,-1,-1,1)^{ \mathrm{\scriptscriptstyle T} }$, and $\psi=1$. This setup satisfies the necessary identifying assumptions, and since $Y^1-Y^0 = \psi T$ it follows that the LIV curve is linear in the latent threshold $T$ (regardless of the conditioning set $\mathbf{V}$). It can also be shown that the above setup implies
\begin{equation*}
\begin{gathered}
\mathbb{E}(A \mid \mathbf{X},Z) = \Phi\left(\frac{Z-\mathbf{o}ldsymbol\beta^T \mathbf{X}}{\sqrt{2}}\right), \\

\mathbb{E}(Y \mid \mathbf{X},Z) = \psi \left\{ (\mathbf{o}ldsymbol\beta^T \mathbf{X}) \Phi\left(\frac{Z-\mathbf{o}ldsymbol\beta^T \mathbf{X}}{\sqrt{2}}\right) - \sqrt{2} \phi\left(\frac{Z-\mathbf{o}ldsymbol\beta^T \mathbf{X}}{\sqrt{2}}\right) \right\},
\end{gathered}
\end{equation*}
so that the treatment regression $\lambda$ follows a probit model; the outcome regression $\mu$ is more complicated, but $\mathbb{E}(Y \mid \mathbf{X}, Z,A)$ follows a particular generalized additive model. 

We fit the proposed estimator, which in this case has a closed form given by
$$ \hat\psi = \frac{ \mathbb{P}_n \left[ \{ w(Z) + Zw'(Z) \} \left\{ \frac{Y-\hat\mu(\mathbf{X},Z)}{\hat\pi(Z \mid \mathbf{X})} \right\} \right]  + \int_\mathcal{T} \{w(t) + tw'(t)\} \hat{m}(t) \ dt}
{ \mathbb{P}_n \left[ \{ 2Z w(Z) + Z^2 w'(Z) \} \left\{ \frac{A-\hat\lambda(\mathbf{X},Z)}{\hat\pi(Z \mid \mathbf{X})} \right\} \right]  + \int_\mathcal{T} \{ 2t w(t) + t^2 w'(t) \} \hat\ell(t) \ dt } . $$
In particular we considered inverse-probability (IP) -weighted estimators with $\hat\mu=\hat\lambda=0$, regression estimators with $\hat\pi=\infty$, and doubly robust estimators that rely on estimates of all of $(\pi,\lambda,\mu)$. {Note that although we recommend the doubly robust estimators in practice, even these IP and regression estimators are novel and have not been previously proposed, to the best of our knowledge}. For the weight function $w(t)$ we used the density of a $\text{TN}\{0,1,(-1.9,1.9)\}$ truncated normal variable, which roughly matches the marginal distribution of $Z$. To misspecify models for $(\hat\pi,\hat\lambda,\hat\mu)$ we used the \textcite{kang2007demystifying} covariate transformations; to misspecify $\hat\lambda$ we additionally transformed $Z$ to $e^Z$ and used a logit rather than probit model. We used maximum likelihood to estimate $\hat\pi$ and $\hat\lambda$, and a generalized additive model to estimate $\hat\mu$. {Coverage was assessed based on bootstrap confidence intervals, using 100 bootstrap samples.} Results (including bias, standard errors, root mean squared error (RMSE), and coverage) are shown in Figure~\ref{fig:simres}.

\begin{figure}
\caption{Simulation results across 500 simulations: bias (with standard error), RMSE (scaled by $\sqrt{n}
\label{fig:simres}
\end{figure}

The weighting estimator gives large bias unless its nuisance estimator $\hat\pi$ is correctly specified, and similarly the regression estimator gives large bias unless $\hat\lambda$ and $\hat\mu$ are correctly specified; however, the doubly robust estimator gives small bias as long as either $\hat\pi$ or $(\hat\lambda,\hat\mu)$ are correct. In our setup, even under misspecification, the bias for the doubly robust estimator is smaller than that of the weighting or regression estimators. The weighted estimator is least precise, while the regression estimator is most precise, and can outperform the doubly robust estimator in terms of mean squared error when both $\hat\lambda$ and $\hat\mu$ are estimated well, particularly for smaller sample sizes (indicating a bias-variance trade-off). {Coverage roughly coincided with bias: the weighting and regression estimators gave poor coverage unless their nuisance estimators were correctly specified, but the doubly robust estimator gave good coverage as long as either $\hat\pi$ or $(\hat\lambda,\hat\mu)$ was modeled well. Correctly-specified weighting and doubly robust estimators exhibited some conservative behavior (i.e., greater than 95\% coverage) at $n=1000$, and slightly anti-conservative behavior for the larger sample size; but this may be explained by the relatively modest number of bootstrap samples and simulations (which was required to reduce computation time).}

{The regression estimator gave slight bias and relatively poor coverage even for large sample sizes; we hypothesize that this is because this estimator relies on a generalized additive model estimator, which may not be smooth enough to guarantee $\sqrt{n}$ rates and asymptotic linearity. As noted throughout our paper, this is an important motivation for using doubly robust and other influence-function-based estimators, which can attain $\sqrt{n}$ rates and asymptotic linearity even if their nuisance estimators do not. In the Supplementary Materials we present simulation results based on using flexible methods (generalized additive models and random forests) to estimate all of $(\pi,\lambda,\mu)$, not just $\mu$. }

\section{Illustration}

In this section we apply the proposed methodology to estimate the effects on infant mortality of delivery at hospitals with high- versus low-level neonatal intensive care units (NICUs). Following \textcite{lorch2012differential} and others, we define high-level NICUs as those that are designated as level III by the American Academy of Pediatrics, and that deliver at least 50 low birthweight infants on average per year. Level III units have high technical capacity, providing subspecialist teams, advanced imaging, and the ability for sustained mechanical assisted ventilation. On the other hand, level I-II NICUs are only designed to provide basic care to lower-risk infants. The question of whether and how care at high-level NICUs might impact infant mortality is important from both patient and policy perspectives. For example if high-level units can reduce infant mortality, particularly among high-risk infants, then policies that send high-risk infants to high-level NICUs might be worth pursuing. 

To assess potential benefits of delivery at hospitals with high-level units, \textcite{lorch2012differential} collected data on all $n=192,078$ premature births in Pennsylvania between 1995 and 2006. Covariate information included data about the infant, such as birthweight and gestational age, as well as about the delivering mother, such as age, race, and measures of socioeconomic status and comorbidities. {A full list of covariate information is given in the Supplementary Materials, and} more details can be found in \textcite{baiocchi2010building} and \textcite{lorch2012differential}. Importantly, the data are missing some detailed clinical information (e.g., comorbidity severity and lab results) that might explain mothers' deliveries at high- versus low-level hospitals; therefore analyses relying on `no unmeasured confounding' assumptions could be suspect. Fortunately,  \textcite{baiocchi2010building} and \textcite{lorch2012differential} identified a potential IV, which is the excess travel time (in minutes) it takes a mother to get to the nearest high- versus low-level NICU. This is a plausible instrument since it affects where mothers deliver (larger values mean mothers have to travel longer to get to high-level units), but it likely does not independently affect infant mortality and is probably not associated with unmeasured confounders that also affect mortality (at least conditional on measured factors like socioeconomic status). More discussion can be found in \textcite{baiocchi2010building} and \textcite{lorch2012differential}. Figure \ref{fig:data} shows loess fits of the unadjusted relationship between instrument and treatment (which is strong), and between  instrument and outcome (which is less strong); the raw data points also indicate the marginal distribution of the instrument. The gray regions denote pointwise 95\% confidence intervals. 

\begin{figure}
\caption{Relationship between instrument $Z$ (excess travel time) and treatment $A$ (delivery at low-level unit) on the left, and instrument and outcome $Y$ (infant mortality) on the right. Note that $y$-axis scaling differs in the two plots. \label{fig:data}
\label{fig:data}
\end{figure}

We conducted two sets of analyses based on the methodology proposed in previous sections. First we estimated the LIV curve only conditional on the threshold value (so that $\mathbf{V}=\emptyset$), and used the proposed cross-validation approach to select among spline models. Second we estimated how effects vary with birthweight and gestational age, which are two important potential effect modifiers. In both analyses it is first necessary to estimate the nuisance functions, which we did using generalized additive models. To estimate the instrument density $\pi$, we used a model previously used by \textcite{kennedy2017nonparametric}, in which the density only depends on covariates through the mean and variance functions but is otherwise flexible. Specifically this model assumes $Z = \pi_1(\mathbf{X}) + \pi_2(\mathbf{X})\epsilon$, where $\epsilon$ satisfies $\mathbb{E}(\epsilon \mid \mathbf{X})=0$ and $\mathbb{E}(\epsilon^2 \mid \mathbf{X})=1$, the density $f_\epsilon$ of $\epsilon$ is unspecified but smooth, and $(\pi_1,\pi_2)$ follow generalized additive models with identity and log links, respectively. Thus under this model the conditional density of the instrument is given by $\pi(z \mid \mathbf{x}) = f_\epsilon[\{ z-\pi_1(\mathbf{x})\} /\pi_2(\mathbf{x}) ]$. 

In our first analysis we estimated the LIV curve $\gamma(t)$ using a density-weighted projection based on the marginal density of the instrument, so that $w(t)=\hat{p}(Z=t)$ for $\hat{p}$ a usual kernel density estimator. {A plot of this weight function is given in the Supplementary Materials. The weighted average treatment effects we report should therefore be interpreted as averages of the true LIV curve values, where the averaging weights non-extreme moderate instrument values most highly, according to the weight function plotted in the Supplementary Materials. As noted in Section 3.2, when projecting onto a constant, the resulting weighted average effect also equals the complier average effect in a population where the density of the latent threshold $T$ equals the weight function $w^*$.} We used natural cubic splines for $\gamma_k(t;\boldsymbol\psi_k)$ with degrees of freedom $k \in \{1,2,3,4\}$ selected via cross-validation with two folds, using the proposed doubly robust pseudo-risk estimate $\hat{R}(\hat\gamma_k)$. 

The linear model with $k=2$ gave the smallest pseudo-risk ($-13.7 \times 10^{-6}$), much smaller than $k=3$ and $k=4$, which led to overfitting ($-7.3 \times 10^{-6}$ and $1.5 \times 10^{-3}$). The risk under $k=2$ was similar to that of the constant effect model with $k=1$ ($-12.6 \times 10^{-6}$) and gave very similar estimates. For example, for the linear model the effect estimates range from 9.0 to 8.9 deaths per 1000 births for excess travel times ranging from 0 to 100, and at level 0.05 we cannot reject the hypothesis that the slope parameter equals zero ($p=0.98$). Table \ref{tab:restab1} gives estimates and 95\% confidence intervals (based on the bootstrap) for three estimators using the constant effect working model $\gamma(t;\boldsymbol\psi)=\psi$;  the inverse-probability-weighted estimator only relies on estimating the conditional instrument density $\pi$ (i.e., it plugs in sample averages of $A$ and $Y$ for $\hat\lambda$ and $\hat\mu$), the regression-based estimator only relies on estimating the treatment and outcome regressions $(\lambda,\mu)$ (i.e., it plugs in $\infty$ for $\hat\pi$), and the doubly robust estimator is the proposed approach detailed in Section 3.

\begin{table}[h!]
\caption{Risk difference estimates using weighted average of LIV curve (in terms of deaths prevented per 1000 births). \label{tab:restab1}}
\begin{center}

\begin{tabular}{l|r}
Method & Est (95\% CI) \\
\hline
Inverse-probability-weighted & -4.8 (-17.2, 7.6) \\
Regression-based & 9.2 (6.3, 12.1)  \\
Doubly robust & 8.9 (5.4, 12.5)  \\
\end{tabular}
\end{center}
\end{table}

The proposed doubly robust estimator indicates a mortality benefit (risk difference) of 8.9 fewer deaths per 1000 births due to high-level NICU care (95\% CI: 5.4, 12.5), among compliers who could be encouraged by travel time to go to a low-level unit. For comparison, this estimate contrasts with the unadjusted risk difference of -18.6 (-20.0, -17.2), which makes high-level NICUs appear to be harming infants, and a doubly robust no-unmeasured-confounding-based estimate of -0.6 (-2.8, 1.6) for the average treatment effect, which does not give any evidence of benefit. Our estimator yields similar results as the two-stage least squares and matching analyses in \textcite{baiocchi2010building}, but targets a different parameter and relies on weaker assumptions. Our regression and doubly robust estimates were similar, and differed somewhat from weighting, indicating that the conditional density model might be misspecified (although the weighted estimator is also imprecise). 

In our second analysis (exploring effect modification by birthweight and gestational age), we projected onto a model in which effects do not vary with the latent threshold (based on the results of our first analysis) but can vary with normal versus low birthweight (2000+ grams versus $<$2000 grams) and early versus very early gestational age (35--37 weeks versus $\leq$34 weeks). Therefore in this analysis we set $\gamma(t,v;\boldsymbol\psi) = \sum_j \psi_j \mathbbold{1}(v=j)$ where $j \in \{1,2,3,4\}$ indexes the four groups. Results are given in Table \ref{tab:restab2}. 

\begin{table}[h!]
\caption{Effect estimates (95\% confidence intervals) by birthweight and gestational age. \label{tab:restab2}}
\begin{center}

\begin{tabular}{l|rr}
 & \multicolumn{2}{c}{{Gestational age}} \\
{Birthweight} & $\leq 34$ wks & 35--37 wks \\
\hline
$<$2000 g & 58.5 (52.7, 64.3) & 6.1 (2.2, 10.0) \\
2000+ g & 10.1 (2.7, 17.4) & 2.4 (-0.8, 5.6) \\
\end{tabular}
\end{center}
\end{table}

The largest effect of high-level NICU care was for the highest-risk infants with low birthweight and very early gestational age; in particular, for this group, care at high-level NICUs was estimated to yield 58.5 fewer deaths per 1000 births (95\% CI: 52.7, 64.3). Effects in the other two higher-risk groups were relatively similar, with an estimated 6.1 and 10.1 fewer deaths per 1000 births (and both statistically significantly different from zero). For the lower-risk group with higher birthweight and gestational age, care at high-level units was estimated to yield 2.4 fewer deaths per 1000 births (95\% CI: -0.8, 5.6), and we cannot reject a null hypothesis of zero average effect. Results were similar but less pronounced when we used less extreme cutoffs for birthweight (2500+ grams versus $<$2500 grams) and gestational age ($\leq$ 35 weeks versus 36--37 weeks).

\section{Discussion}

In this paper we developed novel semiparametric theory and estimation procedures for a marginal version of the LIV curve, which represents the effect among local compliers who would be encouraged to take treatment at a given threshold value of the instrument but not below. Importantly, in contrast to available methods for estimating the fully conditional LIV curve, our methods have the following advantages: they do not require parametric assumptions (but can still yield parametric root-n rates of convergence), incorporate information about the instrument mechanism, are doubly robust (i.e., still yield consistent estimates under misspecification of either the instrument or treatment/outcome processes), and allow for estimating varying amounts of effect modification. We described the asymptotic properties of our methods under weak empirical process conditions, and also proposed a doubly robust cross-validation approach for model selection. Finally we used the proposed methods to study the effects of care at high-level NICUs on infant mortality, including how such effects are modified by infants' birthweight and gestational age. 

There are a number of direct opportunities for future work based on this research{, beyond those areas mentioned in the main text (e.g., extensions to deal with positivity violations, loss functions beyond $L_2$, further development and assessment of the sample-splitting estimator, nonparametric estimation of the local IV curve via kernel-smoothed projections, post-selection inference, etc.)}. First, it will be very useful to develop computationally efficient software for implementing the proposed methods for general non-linear working models. The methods are computationally demanding due to the need to calculate multiple derivatives and integrals, especially in cases involving complex effect modification. Second, it might be of interest to determine the efficient choice of the weight function $w(t,\mathbf{v})$ for the case where the working model $\gamma(t,\mathbf{v};\boldsymbol\psi)$ is believed to be the true model. A third area of future work is in the application studying the effects of high-level NICU care, where it would be useful to implement more flexible covariate adjustment (e.g., Super Learner) and more complex models (e.g., exploring threshold effects and continuous effect modification). 

\stepcounter{section}
\printbibliography[title={\thesection \ \ \ References}]

\pagebreak

\setcounter{page}{1}

\vspace*{.2in}
\begin{center}
\LARGE \textbf{Supplementary Materials for ``Robust causal inference with continuous  instruments  using the local instrumental variable curve''}
\end{center}

\vspace*{.5in}

\normalsize

\setcounter{section}{0}
\section{Proof of Theorem 1}

First note that  
$$ Z \perp\!\!\!\perp (A^z, Y^z) \mid \mathbf{X}  \implies  Z \perp\!\!\!\perp (A^z, Y^{zA^z}) \mid \mathbf{X} \implies  Z \perp\!\!\!\perp (A^z, Y^{A^z}) \mid \mathbf{X} , $$
where the first implication follows from Assumption 2 (consistency) and the second by Assumption 5 (exclusion restriction). \\

Therefore
\begin{align*}
\mathbb{E}(Y \mid \mathbf{X},Z=z ) &= \mathbb{E}(Y^{A^z} \mid \mathbf{X},Z=z) = \mathbb{E}(Y^{A^z} \mid \mathbf{X}) = \mathbb{E}\{A^z (Y^1 -Y^0) \mid \mathbf{X} \} + \mathbb{E}(Y^0 \mid \mathbf{X}) 
\end{align*}
where the first equality follows from Assumption 2 (consistency), the second since Assumption 4 (unconfoundedness of $Z$) implies $Z \perp\!\!\!\perp (A^z, Y^{A^z}) \mid \mathbf{X}$ under Assumptions 2 (consistency) and 5 (exclusion restriction) as shown above, and the third by rearranging. By the same logic we have
$$ \mathbb{E}(A \mid \mathbf{X},Z=z ) = \mathbb{E}(A^z \mid \mathbf{X},Z=z) = \mathbb{E}(A^z \mid \mathbf{X} ) . $$
Assumption 3 (positivity) allows us to write conditional expectations given $\mathbf{X}$ and $Z$. \\

Therefore, combining the above results gives
$$ \mathbb{E}(Y \mid \mathbf{X},Z=z+\delta )  - \mathbb{E}(Y \mid \mathbf{X},Z=z )  = \mathbb{E}\{(A^{z+\delta}-A^z) (Y^1 -Y^0) \mid \mathbf{X}\} $$
and
$$ \mathbb{E}(A \mid \mathbf{X},Z=z+\delta )  - \mathbb{E}(A \mid \mathbf{X},Z=z )  = \mathbb{E}(A^{z+\delta}-A^z \mid \mathbf{X}) , $$
so that 
\begin{align*}
\mathbb{E}\{ \mathbb{E}(Y \mid \mathbf{X}&,Z=z+\delta )  - \mathbb{E}(Y \mid \mathbf{X},Z=z ) \mid \mathbf{V}\} = \mathbb{E}\{(A^{z+\delta}-A^z) (Y^1 -Y^0) \mid \mathbf{V}\} \\
&= \mathbb{E}(Y^1 -Y^0 \mid \mathbf{V} , A^{z+\delta}>A^z ) P(A^{z+\delta}>A^z \mid \mathbf{V}) \\
&=  \mathbb{E}(Y^1 -Y^0 \mid \mathbf{V} , z < T \leq z+\delta) P(z < T \leq z+\delta \mid \mathbf{V}) 
\end{align*}
and similarly
\begin{align*}
\mathbb{E}\{ \mathbb{E}(A \mid \mathbf{X}&,Z=z+\delta )  - \mathbb{E}(A \mid \mathbf{X},Z=z ) \mid \mathbf{V}\} = \mathbb{E}(A^{z+\delta}-A^z \mid \mathbf{V}\} \\
&= P(A^{z+\delta}>A^z \mid \mathbf{V}) =  P(z < T \leq z+\delta \mid \mathbf{V}) ,
\end{align*}
where the first equalities follow by iterated expectation, the second by Assumption 1 (monotonicity), which implies $A^{z+\delta}-A^z = \mathbbold{1}(A^{z+\delta}>A^z)$, and the third by definition of the latent threshold $T$, i.e., $\{A^{z+\delta}>A^z\} \iff \{ A^{z+\delta}=1, A^z=0 \} \iff \{ T \leq z+\delta, T>z \}$.  \\

Therefore, letting $\gamma(\mathbf{v},t)=\mathbb{E}(Y^1-Y^0 \mid T=t, \mathbf{V}=\mathbf{v})$, we have
\begin{align*}
\lim_{\delta \rightarrow 0}  \frac{1}{\delta} &\mathbb{E}\{ \mathbb{E}(Y \mid \mathbf{X},Z=t+\delta )  - \mathbb{E}(Y \mid \mathbf{X},Z=t ) \mid \mathbf{V}\} \\
&= \lim_{\delta \rightarrow 0} \frac{1}{\delta} \mathbb{E}(Y^1 -Y^0 \mid \mathbf{V} , t \leq T \leq t+\delta) P(t \leq T \leq t+\delta \mid \mathbf{V})  \\
&= \gamma(t,\mathbf{V}) \ \lim_{\delta\rightarrow 0} \frac{1}{\delta} \Big\{ P(T \leq t+\delta \mid  \mathbf{V}) - P(T \leq t \mid \mathbf{V}) \Big\} \\
&= \gamma(t,\mathbf{V}) \ p(T =t \mid \mathbf{V})
\end{align*}
and similarly
\begin{align*}
\lim_{\delta \rightarrow 0} \frac{1}{\delta} &\mathbb{E}\{ \mathbb{E}(A \mid \mathbf{X},Z=t+\delta )  - \mathbb{E}(A \mid \mathbf{X},Z=t ) \mid \mathbf{V}\} = \lim_{\delta \rightarrow 0} \frac{1}{\delta} P(t \leq T \leq t+\delta \mid \mathbf{V})  \\
&= \lim_{\delta\rightarrow 0} \frac{1}{\delta} \Big\{ P(T \leq t+\delta \mid \mathbf{V}) - P(T \leq t \mid \mathbf{V}) \Big\} \\
&= p(T =t \mid \mathbf{V}) ,
\end{align*}
where the equalities follow by Assumption 7 (continuity). Specifically, the first and third equalities follow by the fact that $T$ is continuously distributed, with $p(T=t \mid \mathbf{V}) = \frac{\partial}{\partial t} P(T \leq t \mid \mathbf{V})$, and the second follows by the continuity of $\gamma(\mathbf{v},t)$ in $t$. \\

Therefore
$$ \gamma(t,\mathbf{v}) =  \frac{ \frac{\partial}{\partial z} \mathbb{E}\{ \mathbb{E}(Y \mid \mathbf{X},Z=z ) \mid \mathbf{V}=\mathbf{v}\}  }{ \frac{\partial}{\partial z} \mathbb{E}\{ \mathbb{E}(A \mid \mathbf{X},Z=z ) \mid \mathbf{V}=\mathbf{v}\} } \bigg|_{z=t} $$

since the denominator is bounded away from zero by Assumption 6 (instrumentation). \\

\section{Proof of Theorem 2}

In this section we use subscripts to index quantities that depend on the distribution $P$; a zero subscript denotes a quantity evaluated at the true distribution $P=P_0$. Thus for example $\mathbb{E}_P$ denotes expectations under $P$ and $\mathbb{E}_0$ denotes expectations under the truth $P=P_0$; similarly $\boldsymbol\psi_P$ denotes the parameter $\boldsymbol\psi=\boldsymbol\psi(P)$ as a map $\boldsymbol\psi: P \mapsto \mathbb{R}^q$ and $\boldsymbol\psi_0$ denotes its true value evaluated at $P_0$. \\

In the interest of notational simplicity we make two slight abuses of notation. First, in the main text $\mathbf{o}ldsymbol\mathbb{V}phi$ was only proportional to the efficient influence function (i.e., $\mathbf{o}ldsymbol\mathbb{V}phi$ did not include the constant matrix scaling factor, which is unnecesary for solving estimating equations), whereas in this appendix we use $\mathbf{o}ldsymbol\mathbb{V}phi$ to denote the full efficient influence function (including the constant matrix scaling factor). Second, the functions $\mathbf{g}_j(t,\mathbf{v};\boldsymbol\psi)$ (for $j=1,2$) in the main text are denoted by $\mathbf{g}'_j(t,\mathbf{v};\boldsymbol\psi)$ in this appendix. \\

We will show that $\mathbf{o}ldsymbol\mathbb{V}phi(\mathbf{O};\boldsymbol\psi_P,\mathbf{o}ldsymbol\eta_P)=\mathbf{o}ldsymbol\mathbb{V}phi_P(\mathbf{O})$ is the efficient influence function by showing that it is the canonical gradient of the pathwise derivative of $\boldsymbol\psi_P$, i.e., that $\mathbf{o}ldsymbol\mathbb{V}phi_P$ satisfies
$$ \frac{\partial \boldsymbol\psi_\epsilon}{\partial \epsilon} \Bigm|_{\epsilon=0}  = \mathbb{E}_0\{ \mathbf{o}ldsymbol\mathbb{V}phi_0(\mathbf{O}) s_0(\mathbf{O}) \}  $$
where $\boldsymbol\psi_\epsilon = \boldsymbol\psi(P_\epsilon)$ denotes the parameter $\boldsymbol\psi$ evaluated at any regular parametric submodel $\{P_\epsilon : \epsilon\}$ passing through $P_0$ at $\epsilon=0$, and $s_\epsilon(\mathbf{o}_1 \mid \mathbf{o}_2) = \frac{\partial}{\partial\epsilon^*} \log dP_\epsilon^*(\mathbf{o}_1 \mid \mathbf{o}_2)  |_{\epsilon^*=\epsilon}$ denotes the parametric submodel score for any partition $(\mathbf{O}_1, \mathbf{O}_2) \subseteq \mathbf{O}$. \\

By definition we have
\begin{align*}
\boldsymbol\psi_P &= \argmin_{\boldsymbol\psi \in \mathbb{R}^q} \int_{\mathcal{V}} \int_{\mathcal{T}} w(t,\mathbf{v}) \Big\{ \gamma_P(t,\mathbf{v})-\gamma(t,\mathbf{v};\boldsymbol\psi) \Big\}^2 p(T=t \mid \mathbf{v}) \ dt \ dP(\mathbf{v}) 
\end{align*}
and thus
$$ \int_{\mathcal{V}} \int_{\mathcal{T}} \frac{\partial \gamma(t,\mathbf{v};\boldsymbol\psi)}{\partial \boldsymbol\psi}\Bigm|_{\boldsymbol\psi=\boldsymbol\psi_P} w(t,\mathbf{v}) \Big\{ \gamma_P(t,\mathbf{v})-\gamma(t,\mathbf{v};\boldsymbol\psi_P) \Big\} p(T=t \mid \mathbf{v}) \ dt \ dP(\mathbf{v}) = 0 . $$ 

Letting $m_P(z,\mathbf{v})=\mathbb{E}_P\{\mathbb{E}_P(Y \mid \mathbf{X}, Z=z) \mid \mathbf{V}=\mathbf{v}\}$ and $m_P'(t,\mathbf{v})=\frac{\partial}{\partial z} m_P(z,\mathbf{v})|_{z=t}$, and similarly $\ell_P(z,\mathbf{v})=\mathbb{E}_P\{\mathbb{E}_P(A \mid \mathbf{X}, Z=z) \mid \mathbf{V}=\mathbf{v}\}$ and $\ell_P'(t,\mathbf{v})=\frac{\partial}{\partial z} \ell_P(z,\mathbf{v})|_{z=t}$, then under the identifying assumptions in the main text we have
$$ \gamma_P(t,\mathbf{v}) = \frac{m_P'(t,\mathbf{v})}{\ell_P'(t,\mathbf{v})} \ \text{ and } \ p(t \mid \mathbf{v}) = \ell_P'(t,\mathbf{v}) . $$

Therefore the restriction above is equivalent to
\begin{align*}
0 &= \int_{\mathcal{V}} \int_{\mathcal{T}} \frac{\partial \gamma(t,\mathbf{v};\boldsymbol\psi)}{\partial \boldsymbol\psi}\Bigm|_{\boldsymbol\psi=\boldsymbol\psi_P} w(t,\mathbf{v}) \Big\{ m_P'(t,\mathbf{v})  - \gamma(t,\mathbf{v};\boldsymbol\psi_P) \ell_P'(t,\mathbf{v}) \Big\} \ dt \ dP(\mathbf{v}) \\
&= \int_{\mathcal{V}} \int_{\mathcal{T}} \Big\{ \mathbf{g}_2(t,\mathbf{v};\boldsymbol\psi_P) \ m_P'(t,\mathbf{v})  - \mathbf{g}_1(t,\mathbf{v};\boldsymbol\psi_P) \ \ell_P'(t,\mathbf{v}) \Big\} \ dt \ dP(\mathbf{v})
\end{align*}
where $\mathbf{g}_1$ and $\mathbf{g}_2$ are $q$-vectors (with known functional form not depending on $P$) defined as
$$ \mathbf{g}_1(t,\mathbf{v};\boldsymbol\psi) = \mathbf{g}_2(t,\mathbf{v};\boldsymbol\psi) \gamma(t,\mathbf{v};\boldsymbol\psi) \ \text{ and }  \ \mathbf{g}_2(t,\mathbf{v};\boldsymbol\psi) =  \frac{\partial \gamma(t,\mathbf{v};\boldsymbol\psi^*)}{\partial \boldsymbol\psi^*}\Bigm|_{\boldsymbol\psi^*=\boldsymbol\psi} w(t,\mathbf{v})  . $$

And since the weight satisfies $w(t,\mathbf{v})=0$ for $t \notin \text{int}(\mathcal{T})$,  integration by parts gives
$$ \int_{\mathcal{V}} \int_{\mathcal{T}} \Big\{ \mathbf{g}'_1(t,\mathbf{v};\boldsymbol\psi_P) \ \ell_P(t,\mathbf{v})  - \mathbf{g}'_2(t,\mathbf{v};\boldsymbol\psi_P) \ m_P(t,\mathbf{v})  \Big\} \ dt \ dP(\mathbf{v}) = 0 , $$
where $\mathbf{g}'_j(t,\mathbf{v};\boldsymbol\psi)=\partial \mathbf{g}_j(z,\mathbf{v};\boldsymbol\psi)/\partial z |_{z=t}$. \\

Evaluating the above at $P=P_\epsilon$ gives 
$$ \int_{\mathcal{V}} \int_{\mathcal{T}} \Big\{ \mathbf{g}'_1(t,\mathbf{v};\boldsymbol\psi_\epsilon) \ell_\epsilon(t,\mathbf{v})  - \mathbf{g}'_2(t,\mathbf{v};\boldsymbol\psi_\epsilon) m_\epsilon(t,\mathbf{v})  \Big\} \ dt \ dP_\epsilon(\mathbf{v}) = 0 , $$
and differentiating with respect to $\epsilon$ and evaluating at the truth $\epsilon=0$ (using the chain rule) gives 
\begin{align*}
0 &= \int_{\mathcal{V}} \int_{\mathcal{T}} \bigg\{ \frac{\partial \mathbf{g}'_1(t,\mathbf{v};\boldsymbol\psi)}{\partial \boldsymbol\psi} \Bigm|_{\boldsymbol\psi=\boldsymbol\psi_0} \frac{\partial \boldsymbol\psi_\epsilon}{\partial \epsilon}\Bigm|_{\epsilon=0} \ell_0(t,\mathbf{v}) + \mathbf{g}'_1(t,\mathbf{v};\boldsymbol\psi_0) \frac{\partial \ell_\epsilon(t,\mathbf{v})}{\partial \epsilon}\Bigm|_{\epsilon=0} \\
& \hspace{.5in} - \frac{\partial \mathbf{g}'_2(t,\mathbf{v};\boldsymbol\psi)}{\partial \boldsymbol\psi} \Bigm|_{\boldsymbol\psi=\boldsymbol\psi_0} \frac{\partial \boldsymbol\psi_\epsilon}{\partial \epsilon}\Bigm|_{\epsilon=0} m_0(t,\mathbf{v}) - \mathbf{g}'_2(t,\mathbf{v};\boldsymbol\psi_0) \frac{\partial m_\epsilon(t,\mathbf{v})}{\partial \epsilon}\Bigm|_{\epsilon=0}  \bigg\} dt \ dP_0(\mathbf{v}) \\
& \hspace{.25in} + \int_{\mathcal{V}} \int_{\mathcal{T}} \Big\{ \mathbf{g}'_1(t,\mathbf{v};\boldsymbol\psi_0) \ell_0(t,\mathbf{v})  - \mathbf{g}'_2(t,\mathbf{v};\boldsymbol\psi_0) m_0(t,\mathbf{v})  \Big\} s_0(\mathbf{v}) \ dt \ dP_0(\mathbf{v}) .
\end{align*}

Rearranging, this implies that
\begin{align*}
\frac{\partial \boldsymbol\psi_\epsilon}{\partial \epsilon}\Bigm|_{\epsilon=0} &= \mathbf{C}_0^{-1} \int_{\mathcal{V}} \int_{\mathcal{T}} \bigg[ \mathbf{g}'_1(t,\mathbf{v};\boldsymbol\psi_0) \left\{ \frac{\partial \ell_\epsilon(t,\mathbf{v})}{\partial \epsilon}\Bigm|_{\epsilon=0}  + \ell_0(t,\mathbf{v}) s_0(\mathbf{v}) \right\} \\
& \hspace{.5in} - \mathbf{g}'_2(t,\mathbf{v};\boldsymbol\psi_0) \left\{ \frac{\partial m_\epsilon(t,\mathbf{v})}{\partial \epsilon}\Bigm|_{\epsilon=0} + m_0(t,\mathbf{v}) s_0(\mathbf{v}) \right\} \bigg] dt \ dP_0(\mathbf{v})
\end{align*}
with
\begin{align*}
\mathbf{C}_P = & -\int_{\mathcal{V}} \int_{\mathcal{T}} \left\{ \frac{ \partial \mathbf{g}'_1(t,\mathbf{v};\boldsymbol\psi) }{\partial \boldsymbol\psi}\Bigm|_{\boldsymbol\psi=\boldsymbol\psi_P} \ell_P(t,\mathbf{v}) - \frac{ \partial \mathbf{g}'_2(t,\mathbf{v};\boldsymbol\psi) }{\partial \boldsymbol\psi}\Bigm|_{\boldsymbol\psi=\boldsymbol\psi_P} m_P(t,\mathbf{v}) \right\} dt \ dP(\mathbf{v}) ,
\end{align*}
and
\begin{align*}
\frac{\partial \ell_\epsilon(z,\mathbf{v})}{\partial \epsilon}\Bigm|_{\epsilon=0} &= \frac{\partial}{\partial \epsilon} \mathbb{E}_\epsilon \{ \mathbb{E}_\epsilon (A \mid \mathbf{X}, Z=z) \mid \mathbf{V}=\mathbf{v}\} |_{\epsilon=0} = \frac{\partial}{\partial \epsilon} \int_\mathcal{W} \sum_{a \in \{0,1\}} a \ p_\epsilon(a \mid \mathbf{x}, z) \  dP_\epsilon(\mathbf{w} \mid \mathbf{v}) \Big|_{\epsilon=0} \\
&= \int_\mathcal{W} \sum_{a \in \{0,1\}} a \ \Big\{ s_0(a \mid \mathbf{x}, z) + s_0(\mathbf{w} \mid \mathbf{v}) \Big\} \ p_0(a \mid \mathbf{x}, z)  \ dP_0(\mathbf{w} \mid \mathbf{v}) \\
&= \mathbb{E}_0\Big(\mathbb{E}_0\Big[ A \{ s_0(A \mid \mathbf{X}, Z) + s_0(\mathbf{W} \mid \mathbf{V}) \} \Bigm| \mathbf{X}, Z=z \Big] \Bigm| \mathbf{V}=\mathbf{v}\Big) ,
\end{align*}
and by the same logic
$$ \frac{\partial m_\epsilon(z,\mathbf{v})}{\partial \epsilon}\Bigm|_{\epsilon=0} = \mathbb{E}_0\Big(\mathbb{E}_0\Big[ Y \{ s_0(Y \mid \mathbf{X}, Z) + s_0(\mathbf{W} \mid \mathbf{V}) \} \Bigm| \mathbf{X}, Z=z \Big] \Bigm| \mathbf{V}=\mathbf{v}\Big) . $$ 

Now we turn to $\mathbb{E}_0\{ \mathbf{o}ldsymbol\mathbb{V}phi_0(\mathbf{O}) s_0(\mathbf{O}) \}$. The putative efficient influence function $\mathbf{o}ldsymbol\mathbb{V}phi_P$ from the main text is given by
\begin{align*}
\mathbf{o}ldsymbol\mathbb{V}phi_P(\mathbf{O}) &= \mathbf{C}_P^{-1}  \bigg[ \mathbf{g}'_1(Z,\mathbf{V};\boldsymbol\psi_P) \bigg\{ \frac{A - \mathbb{E}_P(A \mid \mathbf{X},Z)}{p(Z \mid \mathbf{X})} \bigg\} - \mathbf{g}'_2(Z,\mathbf{V};\boldsymbol\psi_P) \bigg\{ \frac{Y - \mathbb{E}_P(Y \mid \mathbf{X},Z)}{p(Z \mid \mathbf{X})} \bigg\} \\
& \hspace{.5in} + \int_{\mathcal{T}} \Big\{ \mathbf{g}'_1(t,\mathbf{V};\boldsymbol\psi_P) \mathbb{E}_P(A \mid \mathbf{X},Z=t) - \mathbf{g}'_2(t,\mathbf{V};\boldsymbol\psi_P) \mathbb{E}_P(Y \mid \mathbf{X}, Z=t) \Big\} \ dt \bigg] 
\end{align*}
(note the inclusion of the scaling matrix $\mathbf{C}_P^{-1}$), and $s_0(\mathbf{O})$ is the parametric submodel score, which can be decomposed as
$$ s_0(\mathbf{O}) = s_0(Y, A \mid \mathbf{X}, Z) + s_0(Z \mid \mathbf{X}) + s_0(\mathbf{W} \mid \mathbf{V}) + s_0(\mathbf{V}) . $$ 

Therefore
\begin{align*}
\mathbf{C}_0 &\mathbb{E}_0\Big[ \mathbf{o}ldsymbol\mathbb{V}phi_0(\mathbf{O}) \{ s_0(Y, A \mid \mathbf{X}, Z) + s_0(Z \mid \mathbf{X}) + s_0(\mathbf{W} \mid \mathbf{V}) + s_0(\mathbf{V}) \} \Big] \\
&= \mathbb{E}_0 \bigg[ \mathbf{g}'_1(Z,\mathbf{V};\boldsymbol\psi_0) \bigg\{ \frac{A s_0(A \mid \mathbf{X}, Z) }{p_0(Z \mid \mathbf{X})} \bigg\}  - \mathbf{g}'_2(Z,\mathbf{V};\boldsymbol\psi_0) \bigg\{ \frac{Y s_0(Y \mid \mathbf{X}, Z) }{p_0(Z \mid \mathbf{X})} \bigg\} \\
& \hspace{.5in} + \int_{\mathcal{T}} \Big\{ \mathbf{g}'_1(t,\mathbf{V};\boldsymbol\psi_0) \mathbb{E}_0(A \mid \mathbf{X},Z=t) - \mathbf{g}'_2(t,\mathbf{V};\boldsymbol\psi_0) \mathbb{E}_0(Y \mid \mathbf{X}, Z=t) \Big\} \ dt \\
& \hspace{1in} \times \Big\{ s_0(\mathbf{W} \mid \mathbf{V}) + s_0(\mathbf{V}) \Big\} \bigg] \\
&=\int_\mathcal{V} \int_{\mathcal{T}} \bigg( \mathbf{g}'_1(t,\mathbf{v};\boldsymbol\psi_0) \mathbb{E}_0\Big( \mathbb{E}_0\Big[ A \{ s_0(A \mid \mathbf{X}, Z) + s_0(\mathbf{W} \mid \mathbf{V}) \} \Bigm| \mathbf{X}, Z=t \Big] \Bigm| \mathbf{V}=\mathbf{v} \Big) \\
& \hspace{.5in} - \mathbf{g}'_2(t,\mathbf{v};\boldsymbol\psi_0) \mathbb{E}_0\Big( \mathbb{E}_0\Big[ Y \{ s_0(Y \mid \mathbf{X}, Z) + s_0(\mathbf{W} \mid \mathbf{V}) \} \Bigm| \mathbf{X}, Z=t \Big] \Bigm| \mathbf{V}=\mathbf{v} \Big) \\
& \hspace{.5in} + \Big\{ \mathbf{g}'_1(t,\mathbf{v};\boldsymbol\psi_0) \mathbb{E}_0(A \mid \mathbf{X},Z=t) - \mathbf{g}'_2(t,\mathbf{v};\boldsymbol\psi_0) \mathbb{E}_0(Y \mid \mathbf{X}, Z=t) \Big\} s_0(\mathbf{v}) \bigg) \ dt \ dP_0(\mathbf{v}) \\
&= \int_{\mathcal{V}} \int_{\mathcal{T}} \bigg[ \mathbf{g}'_1(t,\mathbf{v};\boldsymbol\psi_0) \bigg\{ \frac{\partial \ell_\epsilon(z,\mathbf{v})}{\partial \epsilon}\Bigm|_{\epsilon=0}  + \ell_0(t,\mathbf{v}) s_0(\mathbf{v}) \bigg\} \\
& \hspace{.5in} - \mathbf{g}'_2(t,\mathbf{v};\boldsymbol\psi_0) \bigg\{ \frac{\partial m_\epsilon(t,\mathbf{v})}{\partial \epsilon}\Bigm|_{\epsilon=0} + m_0(t,\mathbf{v}) s_0(\mathbf{v}) \bigg\} \bigg] dt \ dP_0(\mathbf{v}) = \mathbf{C}_0 \frac{\partial \boldsymbol\psi_\epsilon}{\partial \epsilon}\Bigm|_{\epsilon=0}
\end{align*}
where the first equality follows by iterated expectation and the fact that $\mathbb{E}_0\{ s_0(\mathbf{O}_1 \mid \mathbf{O}_2) \mid \mathbf{O}_2\}=0$ for any $(\mathbf{O}_1,\mathbf{O}_2) \subseteq \mathbf{O}$, the second follows by iterated expectation, the third follows by iterated expectation and by definition of $\ell_0$ and $m_0$ (along with the earlier results for their derivatives with respect to $\epsilon$), and the fourth follows by the expression derived earlier for ${\partial \boldsymbol\psi_\epsilon} / {\partial \epsilon}|_{\epsilon=0}$. \\

Therefore, as long as $\mathbf{C}_0$ is invertible, we have ${\partial \boldsymbol\psi_\epsilon} / {\partial \epsilon}|_{\epsilon=0} = \mathbb{E}_0\{ \mathbf{o}ldsymbol\mathbb{V}phi_0(\mathbf{O}) s_0(\mathbf{O}) \}$ and thus  $\mathbf{o}ldsymbol\mathbb{V}phi_P(\mathbf{O})$ is the efficient influence function.

\section{Double robustness of efficient influence function $\mathbf{o}ldsymbol\mathbb{V}phi$}

Here we will show that $\mathbb{E}\{ \mathbf{o}ldsymbol\mathbb{V}phi(\mathbf{O}; \boldsymbol\psi,\overline\pi,\overline\lambda,\overline\mu)\}=0$ as long as either 
$$\overline\pi=\pi_0 \ \text{ or } \ (\overline\lambda,\overline\mu)=(\lambda_0,\mu_0). $$ 
In this section expectations $\mathbb{E}=\mathbb{E}_0$ and parameters $\boldsymbol\psi=\boldsymbol\psi_0$ are evaluated under $P_0$, but we drop the subscript for notational convenience. \\

First note that
\begin{align*}
\mathbf{C}_0 & \mathbb{E}\{ \mathbf{o}ldsymbol\mathbb{V}phi(\mathbf{O}; \boldsymbol\psi,\overline\pi,\overline\lambda,\overline\mu)\} = \bigg[ \mathbf{g}'_1(Z,\mathbf{V};\boldsymbol\psi) \bigg\{ \frac{A - \overline\lambda(\mathbf{X},Z)}{\overline\pi(Z \mid \mathbf{X})} \bigg\} - \mathbf{g}'_2(Z,\mathbf{V};\boldsymbol\psi) \bigg\{ \frac{Y - \overline\mu(\mathbf{X},Z)}{\overline\pi(Z \mid \mathbf{X})} \bigg\} \\
& \hspace{1.7in} + \int_{\mathcal{T}} \Big\{ \mathbf{g}'_1(t,\mathbf{V};\boldsymbol\psi) \overline\lambda(\mathbf{X},t) - \mathbf{g}'_2(t,\mathbf{V};\boldsymbol\psi) \overline\mu(\mathbf{X}, t) \Big\} \ dt \bigg] \\
&= \mathbb{E}\bigg[ \mathbf{g}'_1(Z,\mathbf{V};\boldsymbol\psi) \bigg\{ \frac{\lambda_0(\mathbf{X},Z)- \overline\lambda(\mathbf{X},Z)}{\overline\pi(Z \mid \mathbf{X})} \bigg\} - \mathbf{g}'_2(Z,\mathbf{V};\boldsymbol\psi) \bigg\{ \frac{\mu_0(\mathbf{X},Z) - \overline\mu(\mathbf{X},Z)}{\overline\pi(Z \mid \mathbf{X})} \bigg\} \\
& \hspace{.5in} + \int_{\mathcal{T}} \Big\{ \mathbf{g}'_1(t,\mathbf{V};\boldsymbol\psi) \overline\lambda(\mathbf{X},t) - \mathbf{g}'_2(t,\mathbf{V};\boldsymbol\psi) \overline\mu(\mathbf{X}, t) \Big\} \ dt \bigg] \\
&= \mathbb{E}\int_{\mathcal{Z}} \bigg[ \mathbf{g}'_1(t,\mathbf{V};\boldsymbol\psi) \Big\{ \lambda_0(\mathbf{X},t)- \overline\lambda(\mathbf{X},t) \Big\}  - \mathbf{g}'_2(t,\mathbf{V};\boldsymbol\psi) \Big\{ \mu_0(\mathbf{X},t) - \overline\mu(\mathbf{X},t) \Big\} \bigg] \frac{\pi_0(t \mid \mathbf{X})}{\overline\pi(t \mid \mathbf{X})} \ dt \\
& \hspace{.5in} + \int_{\mathcal{T}} \Big\{ \mathbf{g}'_1(t,\mathbf{V};\boldsymbol\psi) \overline\lambda(\mathbf{X},t) - \mathbf{g}'_2(t,\mathbf{V};\boldsymbol\psi) \overline\mu(\mathbf{X}, t) \Big\} \ dt \bigg] \\
&= \mathbb{E}\int_{\mathcal{T}} \bigg[ \mathbf{g}'_1(t,\mathbf{V};\boldsymbol\psi) \Big\{ \lambda_0(\mathbf{X},t)- \overline\lambda(\mathbf{X},t) \Big\}  - \mathbf{g}'_2(t,\mathbf{V};\boldsymbol\psi) \Big\{ \mu_0(\mathbf{X},t) - \overline\mu(\mathbf{X},t) \Big\} \bigg] \bigg\{ \frac{\pi_0(t \mid \mathbf{X})}{\overline\pi(t \mid \mathbf{X})} - 1\bigg\} dt  \\
& \hspace{.5in} + \int_{\mathcal{T}} \Big\{ \mathbf{g}'_1(t,\mathbf{V};\boldsymbol\psi) \lambda_0(\mathbf{X},t) - \mathbf{g}'_2(t,\mathbf{V};\boldsymbol\psi) \mu_0(\mathbf{X}, t) \Big\} \ dt \bigg]
\end{align*}
where the first equality is true by definition, the second and third holds by iterated expectation given $(\mathbf{X},Z)$ and $\mathbf{X}$, respectively, and the last follows after rearranging and since $\mathbf{g}_1'=\mathbf{g}_2'=0$ for $t \notin \text{int}(\mathcal{T})$. \\

Therefore if  $\overline\pi=\pi_0$ or $(\overline\lambda,\overline\mu)=(\lambda_0,\mu_0)$ then $\mathbf{C}_0 \mathbb{E}\{ \mathbf{o}ldsymbol\mathbb{V}phi(\mathbf{O}; \boldsymbol\psi,\overline\pi,\overline\lambda,\overline\mu)\}$ equals
\begin{align*} 
\int_{\mathcal{T}} \mathbb{E}&\Big\{ \mathbf{g}'_1(t,\mathbf{V};\boldsymbol\psi) \lambda_0(\mathbf{X},t) - \mathbf{g}'_2(t,\mathbf{V};\boldsymbol\psi) \mu_0(\mathbf{X}, t) \Big\} \ dt  \\
&= \int_{\mathcal{V}} \int_{\mathcal{T}} \mathbf{g}'_1(t,\mathbf{v};\boldsymbol\psi) \ \ell_0(t),\mathbf{v}) - \mathbf{g}'_2(t,\mathbf{v};\boldsymbol\psi) \ m_0(t,\mathbf{v}) \Big\} \ dt \ dP(\mathbf{v}) = 0
\end{align*}
where the first equality follows by iterated expectation given $\mathbf{V}$ (and by the definitions of $\ell$ and $m$ from the previous section), and the second follows from the restriction given in the previous section (after using integration by parts).

\section{Sample-splitting estimator}

For simplicity, the estimator presented in the main text solves an estimating equation that depends on nuisance estimates $\mathbf{o}ldsymbol{\hat\eta}=(\hat\pi,\hat\lambda,\hat\mu)$ constructed from the entire sample. The asymptotic normality of such estimators requires empirical process (e.g., Donsker) conditions that restrict the complexity of the nuisance functions and their estimators. Unfortunately, these conditions may not be satisfied for very adaptive machine learning methods, such as random forests; however, such conditions are not necessary for estimators that use sample splitting to separate nuisance function estimation from estimating equation evaluation. We present such an estimator here, following \textcite{robins2008higher, zheng2010asymptotic, chernozhukov2016double}. \\

Sample splitting has a long history (cf.\ Bickel (1982), Schick (1986), van der Vaart (1998)), though the idea of explicitly combining with machine learning came more recently. \textcite{robins2008higher} (page 379) used sample splitting to avoid Donsker conditions for functional estimation in nonparametric models, although they did not use cross-fitting. The PhD thesis of Ayygari (2010) (subsequently published as Robins et al.\ (2013)) used cross-fitting in semiparametric efficient estimation with cross-validation-based selection of nuisance estimators. \textcite{zheng2010asymptotic} gave a general cross-fitting version of targeted maximum likelihood estimation to avoid Donsker conditions. Belloni et al.\ (2010) used cross-fitting to relax sparsity requirements in a high-dimensional linear instrumental variable setting. Cross-fitting was also studied more recently by \textcite{chernozhukov2016double}. \\

We rely on the notation presented in Section 3.4 of the main text. Again let $\mathbf{S}=(S_1,...,S_n)$ denote a random variable independent of the sample that splits the data into training ($S_i=0$) and test ($S_i=1$) sets. Here we suppose the split variable $\mathbf{S}$ takes $v$ different values $\{ \mathbf{S}_1,...,\mathbf{S}_v\}$, each with equal probability $1/v$, where $\sum_i S_{iv}=n/v$ for all $v$ and $\sum_v S_{iv}=1$ for all $i$, so that test sets are all of size $n/v$ and each unit is only used in one test set. Define $\mathbb{P}_\mathbf{s}^0$ and $\mathbb{P}_\mathbf{s}^1$ as the sub-empirical distributions for the training data $\{i: S_i=0\}$ and test data $\{i: S_i=1\}$, respectively, for a given split $\mathbf{S}=\mathbf{s}$. Therefore, for example, $\mathbf{o}ldsymbol{\hat\eta}(\mathbb{P}_\mathbf{s}^0)$ denotes the nuisance function estimates based only on the training set data. Then the proposed sample-splitting version of our estimator is given by the solution to
$$ \mathbb{E}_\mathbf{S}  \mathbb{P}_\mathbf{S}^1[ \mathbf{o}ldsymbol\mathbb{V}phi\{\mathbf{O}; \mathbf{o}ldsymbol\psi, \mathbf{o}ldsymbol{\hat\eta}(\mathbb{P}_\mathbf{S}^0) \} ] = \mathbb{E}_\mathbf{S}  \int \mathbf{o}ldsymbol\mathbb{V}phi\{\mathbf{o}; \mathbf{o}ldsymbol\psi, \mathbf{o}ldsymbol{\hat\eta}(\mathbb{P}_\mathbf{S}^0) \} \ d\mathbb{P}_\mathbf{S}^1(\mathbf{o}) = \mathbf{0} ,  $$
with $\mathbf{o}ldsymbol\mathbb{V}phi$ as defined in the main text.

\section{Proof of Theorem 3}

Since $\hat\boldsymbol\psi$ is a Z-estimator, Theorem 3 follows directly from Theorem 5.31 of van der Vaart (2000), together with the fact that $\mathbf{C}_0 \mathbb{P} \{ \mathbf{o}ldsymbol\mathbb{V}phi(\mathbf{O}; \boldsymbol\psi, \hat{\mathbf{o}ldsymbol\eta}) \} = \mathbf{C}_0 \int \mathbf{o}ldsymbol\mathbb{V}phi(\mathbf{o}; \boldsymbol\psi, \hat{\mathbf{o}ldsymbol\eta}) \ dP(\mathbf{o})$ equals 
\begin{align*}
\mathbb{P}\int_{\mathcal{T}} & \bigg[ \mathbf{g}'_1(t,\mathbf{V};\boldsymbol\psi) \Big\{ \lambda_0(\mathbf{X},t)- \hat\lambda(\mathbf{X},t) \Big\}  - \mathbf{g}'_2(t,\mathbf{V};\boldsymbol\psi) \Big\{ \mu_0(\mathbf{X},t) - \hat\mu(\mathbf{X},t) \Big\} \bigg] \bigg\{ \frac{\pi_0(t \mid \mathbf{X})}{\hat\pi(t \mid \mathbf{X})} - 1\bigg\} dt \\
&= \mathbb{P}\left[ \Big\{ \mathbf{g}'_1 ( \lambda_0 - \hat\lambda)  - \mathbf{g}'_2 ( \mu_0 - \hat\mu ) \Big\} ( \pi_0 - \hat\pi ) \Big/ \hat\pi \pi_0 \right] \leq C || \mathbf{g}'_1 ( \lambda_0 - \hat\lambda)  - \mathbf{g}'_2 ( \mu_0 - \hat\mu )  ||_2  || \pi_0 - \hat\pi ||_2 \\
&= O_p\left\{ \Big(  || \lambda_0 - \hat\lambda ||_2 + ||  \mu_0 - \hat\mu  ||_2 \Big) || \pi_0 - \hat\pi ||_2 \right\}
\end{align*}
where the first line follows by iterated expectation (first given $\mathbf{X}$ and $Z$, then given $\mathbf{X}$), the second by multiplying and dividing by $\pi(t \mid \mathbf{X})$, the inequality follows by Cauchy-Schwarz (i.e., $\mathbb{P}(fg) \leq ||f||_2 ||g||_2$) and boundedness of $1/\hat\pi \pi_0$ (from the Donsker condition and positivity), and the last line by the triangle inequality and boundedness of $\mathbf{g}'_1$ and $\mathbf{g}'_2$.

\section{Avoiding Donsker conditions with sample-splitting}

For the asymptotic properties of the sample-splitting estimator presented in Section 4 of the Supplementary Materials, we refer to \textcite{chernozhukov2016double}, who give a thorough analysis of the conditions required for asymptotic linearity.  \\

Conditions (i)--(v) of their Assumption 2.1 are standard regularity and identifiability conditions in Z-estimator problems, which can be verified for sufficiently smooth working models $\gamma(t,\mathbf{v};\boldsymbol\psi)$ and estimating functions with sufficiently well-separated solutions in expectation. Conditions (i)--(vi) of their Assumption 2.2 will be satisfied under the rate conditions given in our Assumption 4 of Theorem 3 in the main text, along with sufficient smoothness of the working models $\gamma(t,\mathbf{v};\boldsymbol\psi)$. Notably, only rate conditions are required for the nuisance estimators $\mathbf{o}ldsymbol{\hat\eta}$, and not Donsker or other regularity conditions. \textcite{zheng2010asymptotic} give similar results for sample-splitting-based targeted maximum likelihood estimators, relying on higher-level conditions.

\section{Proof of Theorem 4 \& double robustness of loss $L$}

After replacing $\mathbf{g}_j$ with $f_j$, these proofs follow precisely the same logic as the proofs given in previous sections of Theorem 2 and of double robustness of the efficient influence function $\mathbf{o}ldsymbol\mathbb{V}phi$, respectively, and so are omitted. 

\section{Additional simulation results}

In the simulations in the main text, we used correct models for the nuisance functions (parametric models for $(\pi,\lambda)$ and a generalized additive model for $\mu$), with misspecification coming from covariate transformations. In this section we present results using generalized additive models and random forests for all nuisance estimators, since it is of interest to explore the performance of our approach when all nuisance functions are estimated without relying on parametric models. The nuisance functions $(\lambda,\mu)$ are just regression functions so we simply use generalized additive models and random forests as usual there; since $\pi$ is a density we first use generalized additive models or random forests to estimate the conditional expectation $\mathbb{E}(Z \mid \mathbf{X})$, and then use a kernel density estimator to estimate the density depending on the sign of the regression function. Note that these approaches are not necessarily correctly specified (except for estimating $\mu$ with a generalized additive model, which actually holds). As in the main text, we compare performance using both the ``true'' covariates $\mathbf{X}$ and the ``misspecified'' covariates based on the Kang \& Schafer (2007) transformations. \\

Results are given in Figure 1 below. Since the models were not necessarily correctly specified, the findings are harder to generalize. The inverse-probability-weighted estimator gave some bias in all settings, which suggests that the density estimation procedure could be improved. The regression estimator tended to give small bias when it relied on the true covariates for the $\mu$ estimator, as might be expected. The doubly robust estimator gave smallest bias when it relied on the true covariates for the $\pi$ estimator; it is unclear why there was bias otherwise. \\

The findings were slightly more consistent for RMSE, especially when using random forests. However, for both generalized additive model and random forests approaches, the doubly robust estimator gave smallest average RMSE across simulation settings, with the inverse-probability-weighted estimator yielding highest average RMSE, and the regression estimator somewhere in between.

\begin{figure}
\caption{Simulation results across 500 simulations: bias (with standard error), RMSE (scaled by $\sqrt{n}
\label{fig:tab}
\end{figure}

\section{Details about empirical example}

\subsection{List of adjustment covariates}

The data analysis in Section 5 of the main text relied on the following 16 adjustment covariates, as discussed in more detail by Baiocchi et al.\ (2010) and Lorch et al.\ (2012). \\

Information about the zip code in which the mother lives (6):
\begin{itemize}
\item median income, percentage below poverty, median home value, percent with high school degree, percent with college degree, percent who rent versus own home.
\end{itemize}

Information about the mother (8):
\begin{itemize}
\item age, diabetes status, month prenatal care was started, number of times previously given birth, whether multiple deliveries, education level (8th grade or less, some high school, high school graduate, some college, college graduate, or more than college), mother's race (White, Black, Asian/Pacific Islander, or other), insurance type (fee for service, HMO, federal/state, other, or uninsured).
\end{itemize}

Information about the infant (2):
\begin{itemize}
\item birthweight, gestational age.
\end{itemize}

\subsection{Plot of estimated density of $Z$ (i.e., weight function $w$)}

\vspace*{-.3in}
\begin{figure}
\caption{Kernel estimate of the marginal density $p(z=t)$ of the instrument $Z$, also used as the weight function $w(t)=\hat{p}
\end{figure}

\end{document}