\tit
Marc Ratkovic\thanks{Assistant Professor
Department of Politics
Princeton University
Princeton NJ 08544. Phone: 608-658-9665
Email: \href{mailto:ratkovic@princeton.edu
The intersection of causal inference and machine learning is a rapidly advancing field. We propose a new approach, the method of direct estimation, that draws on both traditions in order to obtain nonparametric estimates of treatment effects. The approach focuses on estimating the effect of fluctuations in a treatment variable on an outcome. A tensor-spline implementation enables rich interactions between functional bases allowing for the approach to capture treatment/covariate interactions. We show how new innovations in Bayesian sparse modeling readily handle the proposed framework, and then document its performance in simulation and applied examples. Furthermore we show how the method of direct estimation can easily extend to structural estimators commonly used in a variety of disciplines, like instrumental variables, mediation analysis, and sequential $g$-estimation. \bigskip \noindent {Key Words: causal inference, treatment effects, instrumental variables, mediation, controlled direct effects, Bayesian LASSO, sure independence screening
mde_arxiv.tex
\documentclass[12pt,letterpaper]{article}

%% === submission
\newcommand{\blind}{0}

%% === margins ===
\addtolength{\hoffset}{-0.8in} \addtolength{\voffset}{-0.8in}
\addtolength{\textwidth}{1.6in} \addtolength{\textheight}{1.6in}

%% === basic packages ===
\usepackage{latexsym,multirow}
\usepackage{amssymb,amsmath, bm}
\usepackage{graphicx}
%% === bibliography packages ===
\usepackage{natbib}
\bibliographystyle{apsr}
%% === hyperref options ===
\usepackage{color}
\usepackage[pdftex, bookmarksopen=true, bookmarksnumbered=true,
pdfstartview=FitH, breaklinks=true, urlbordercolor={0 1 0}, citebordercolor={0 0 1}]{hyperref}

\usepackage{subfigure}
\usepackage{placeins}
\usepackage{comment}

% === check mark
\usepackage{pifont}
\newcommand{\cmark}{\ding{51}}
\newcommand{\xmark}{\ding{55}}

% === dcolumn package ===
\usepackage{dcolumn}
\newcolumntype{.}{D{.}{.}{-1}}
\newcolumntype{d}[1]{D{.}{.}{#1}}

% === tikz for pictures ===
\usepackage{tikz}
\usepackage[latin1]{inputenc}
\usetikzlibrary{shapes,arrows,trees,fit,positioning}

% === theorem package ===
\usepackage{theorem}
\theoremstyle{plain}
\theoremheaderfont{\scshape}
\newtheorem{proposition}{Proposition}
\newtheorem{assumption}{Assumption}
\newtheorem{lemma}{Lemma}
\newtheorem{definition}{Definition}
\newcommand{\qed}{\hfill \ensuremath{\Box}}
\newcommand{\indep}{\mbox{$\perp\!\!\!\perp$}}
\newcommand{\notindep}{\mbox{$\not\perp\!\!\!\perp$}}
\def\independenT#1#2{\mathrel{\rlap{$#1#2$}\mkern2mu{#1#2}}}
\DeclareMathOperator{\sgn}{sgn}
\DeclareMathOperator{\diag}{diag}
\providecommand{\norm}[1]{\lVert#1\rVert}

% ==== rotating package ===
\usepackage{rotating}

% ==== dotted lines in tables ===
\usepackage{arydshln}

% == spacing between sections and subsections
\usepackage[compact]{titlesec}
%\usepackage{times}


% == TpX ===
%\usepackage{color}
%\usepackage[pdftex]{graphicx}
%\DeclareGraphicsExtensions{.pdf,.png,.mps}
%\usepackage{pgf}
%\usepackage{tikz}
%\usepackage{epic,bez123}
%\usepackage{floatflt}% package for floatingfigure environment
%\usepackage{wrapfig}% package for wrapfigure environment

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}
\fontsize{12pt}{20 pt plus 1 pt minus 0.5 pt}\selectfont

%Special instructions to get rid of the annoying vector notation
\newcommand\vm[1]{% Vector or matrix
  \bm{\mathrm{#1}}}
% === new commands ===
\newcommand\ud{\mathrm{d}}
\newcommand\dist{\buildrel\rm d\over\sim}
\newcommand\ind{\stackrel{\rm indep.}{\sim}}
\newcommand\iid{\stackrel{\rm i.i.d.}{\sim}}
\newcommand\logit{{\rm logit}}
\newcommand\convd{\stackrel{\rm d}{\rightarrow}}
\newcommand\convp{\stackrel{\rm p}{\rightarrow}}
\renewcommand\r{\right}
\renewcommand\l{\left}
\newcommand\pre{{(t-1)}}
\newcommand\cur{{(t)}}
\newcommand\cA{\mathcal{A}}
\newcommand\bone{\mathbf{1}}
\newcommand\E{\mathbb{E}}
\newcommand\Var{{\rm Var}}
\newcommand\Cov{{\rm Cov}}
\newcommand\med{{\rm med}}
\newcommand\cC{\mathcal{C}}
\newcommand\cD{\mathcal{D}}
\newcommand\cJ{\mathcal{J}}
\newcommand\cN{\mathcal{N}}
\newcommand\cP{\mathcal{P}}
\newcommand\cT{\mathcal{T}}
\newcommand\cX{\mathcal{X}}
\newcommand\wX{\widetilde{X}}
\newcommand\wT{\widetilde{T}}
\newcommand\wY{\widetilde{Y}}
\newcommand\wZ{\widetilde{Z}}
\newcommand\oX{\overline{X}}
\newcommand\oT{\overline{T}}
\newcommand\oY{\overline{Y}}
\newcommand\bX{\mathbf{X}}
\newcommand\bB{\mathbf{B}}
\newcommand\bD{\mathbf{D}}
\newcommand\bM{\mathbf{M}}
\newcommand\bH{\mathbf{H}}
\newcommand\bI{\mathbf{I}}
\newcommand\bG{\mathbf{G}}
\newcommand\bR{\mathbf{R}}
\newcommand\bS{\mathbf{S}}
\newcommand\bV{\mathbf{V}}
\newcommand\bW{\mathbf{W}}
\newcommand\plim{\textrm{plim}\;}
\newcommand\supp{\textrm{supp}\;}
\newcommand\card{\textrm{card}\;}

\newcommand{\argmax}{\operatornamewithlimits{argmax}}
\newcommand{\argmin}{\operatornamewithlimits{argmin}}

\newcommand\spacingset[1]{\renewcommand{\baselinestretch}%
{#1}\small\normalsize}

\spacingset{1.25}

\newcommand{\tit}{ Causal Inference through the Method of Direct Estimation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\if0\blind

{\title{\bf \tit \thanks{We would like to thank Horacio Larreguy for providing replication data. We would also like to thank Peter Aronow, James Fowler, Kosuke Imai, Gary King, John Londregan, Chris Lucas, Shiro Kuriwaki, Rich Nielsen, Brandon Stewart, Aaron Strauss, Teppei Yamamoto, and Xiang Zhou, as well as the participants at the Quantitative Social Science Seminar at Princeton  and Yale Research Design and Causal Inference seminar. Not for citation or distribution without permission from the authors.}}

    \author{Marc Ratkovic\thanks{Assistant Professor, Department of
      Politics, Princeton University, Princeton NJ 08544. Phone:
      608-658-9665, Email:
      \href{mailto:ratkovic@princeton.edu}{ratkovic@princeton.edu},
      URL:
      \href{http://www.princeton.edu/~ratkovic}{http://www.princeton.edu/$\sim$ratkovic}}\hspace{.1in}
    Dustin Tingley\thanks{Professor of Government, Harvard University, Email:
      \href{mailto:dtingley@gov.harvard.edu}{dtingley@gov.harvard.edu},
      URL:
      \href{http://scholar.harvard.edu/dtingley}{scholar.harvard.edu/dtingley} } }




\date{\today}

\maketitle

}\fi

\if1\blind
\title{\bf \tit}

}

\maketitle
\fi

\pdfbookmark[1]{Title Page}{Title Page}

\thispagestyle{empty}
\setcounter{page}{0}

\begin{abstract}



The intersection of causal inference and machine learning  is a rapidly advancing field. We propose a new approach, the method of direct estimation, that draws on both traditions in order to obtain nonparametric estimates of treatment effects. The approach focuses on estimating the effect of fluctuations in a treatment variable on an outcome. A tensor-spline implementation enables rich interactions between functional bases allowing for the approach to capture treatment/covariate interactions. We show how new innovations in Bayesian sparse modeling readily handle the proposed framework, and then document its performance in simulation and applied examples. Furthermore we show how the method of direct estimation can easily extend to structural estimators commonly used in a variety of disciplines, like instrumental variables, mediation analysis, and sequential $g$-estimation.




  \bigskip
  \noindent {\bf Key Words:} causal inference, treatment effects, instrumental variables, mediation, controlled direct effects, Bayesian LASSO, sure independence screening
\end{abstract}


\clearpage
\spacingset{1.5}

%% Considered Notes--win on continuous, even in binary case
%% suspect good performance on binary--interactions/splines
%% For any data set--a larger number that are on common support
%% Look at spline bases for support

%%  Aaron--just doing better job at observationalists
%%            Useful for experimentalists/subgroup analysis
%%            Look at spline bases
%%           Overfitting: No!  SIS and regularization through the prior
%%  Peter--econ literature \& engineering literature
%%           Angrist paper--partial effects
%%  Horacio--play up awesomeness of defiers and LICE/relaxing monotonicity
%%           Simulation 1--TSLS returns something, we don't
%%           Multiple instruments
%%  IV sims--More estimated pos than neg (higher values)




\section{Introduction}

We direct our attention to a single question: what is the expected impact of perturbing a treatment variable on an outcome, holding the remaining variables fixed?  This notion of a ceteris paribus manipulation sits at the center of several formulations of causality \citep{holl:86}.  It can be expressed in equivalent terms as a difference in potential outcomes in the Neyman-Rubin model \citep{Imbens2015, Rubin2005}, the result of a ``do'' operation on the treatment in Pearl's framework \citep{pearl1995causal,pearl2009causality}, or the change in the conditional mean function attributable  to the treatment   \citep{Robins94, hanseneconometrics}.



Identification follows from standard assumptions of no omitted confounders, that the manipulation is on the support of the treatment variable, non-interference of units, and a stable treatment value.  Estimation, though, poses a particular problem.  Estimation methods are often bespoke, tailored to a binary treatment, continuous treatment, instrumental variable model, mediation model, and so on.  These methods often only recover an average effect, or, equivalently, low-dimensional structural parameter, averaging over potentially interesting underlying heterogeneities.


Rather than focus on a particular treatment regime, we instead use the observed covariates to directly estimate the outcome and then use the model to predict counterfactual quantities of interest.  We refer to this as the {\it method of direct estimation}, MDE.  As we are predicting unobserved counterfactual values, the method applies equally to a continuous, binary, or categorical treatment regime.   We focus primarily on estimating the effect of a small perturbation of the treatment, a partial derivative, on the outcome.  We show that this estimand is both useful and illuminating, particularly in structural models.  We extend beyond the local perturbation, though.   For example, when we consider the treatment effect on the treated, we fit a single model to the observed outcome and then predict the  outcome for the treated observations had they counterfactually not received the treatment.  Similarly, in an instrumental variable setting, we estimate two covariances: that between the outcome and the instrument and that between the treatment and the instrument, each conditioning on observed covariates.  Our instrumental variables estimate is the observation-level ratio between these two estimates.  Standard instrumental variables estimates are the average of covariances, masking which observations are actually impacted by the instrument.  As we maintain these heterogeneities, we are able to relax the standard assumption that the instrument has a monotonic effect on the treatment over the sample.

The proposed estimation method has four steps.  First, we generate a large number of nonparametric tensor-spline bases, on the order of millions, including all $treatment\times covariate \times covariate$ interactions and lower-order terms.  These bases constitute a larger space than those considered by earlier multivariate nonparametric regression models that inspire our work \citep{Stone1997,friedman1991multivariate}.  We include such a dense parameterization so as to directly model the confounding between treatment and covariate.  We utilize nonparametric bases as we do not know the functional form this interaction may take, and the nature of the interaction may vary from observation to observation.   Second, we use a marginal correlation screen \citep{Fan:2008} to reduce the number of bases from millions to hundreds or thousands.  Third, we regress the outcome onto the maintained bases using a sparse Bayesian regression.  Fourth, we use the estimated model to predict the counterfactual quantities of interest.

Our larger goal is aligning machine learning, particularly high dimensional nonparametric regression, with causal inference, tightly and at a conceptual level.  For example, we show that the ignorability assumption common to causal inference actually implies a local linear function for the outcome, guiding our choice of bases.  Similarly, we show that our proposed regression method satisfies an oracle inequality not only on the fitted values but also on the treatment effect.  This gives us an optimality condition in predicting the local perturbations, which is close but not exactly the same as predicting the outcome well.\footnote{As one example, ensemble tree methods will fit the systematic component well, but since the fits are piecewise constant, they are not well-suited to estimating the effect of a local perturbation.} The oracle bound has been used as a descriptive tool, giving rates of risk convergence.  Instead, we use it constructively, using a feasible version of the bound as a means of estimating which observations are likely not impacted by an instrument in an encouragement design.


%Our statistical advances are threefold. First, we show that our constructed bases produce regression estimates that are comparable to cutting-edge machine learning methods.  Second, we implement a sparse Bayesian modeling designed to generate estimates that satisfy an oracle bound.  The model includes a prior such that the key tuning parameter grows  in the sample size and number of covariates so as to guarantee an oracle bound.  The prior also includes adaptive weights, similar to those in \citet{zouadaptive}, except they are estimated endogenously.  Furthermore, the prior is nonseparable, so a global tuning parameter helps the parameters adapted to the estimated level of overall sparsity in the model \citep[similar to the nonseparable priors in ][]{bhattacharya2015dirichlet, Rockova17}.  In simulations we show that the gains over existing sparse Bayesian models \citep{cps2010horseshoe, polson2014bayesian} are sizable.  The third major advance is theoretical.  Our key tuning parameter is set to grow at an oracle rate.  Estimating an outcome well is not precisely the same as predicting the impact of perturbing the treatment.  As one example, ensemble tree methods will fit the systematic component well, but since the fits are piecewise constant, they are not well-suited to estimating the effect of a local perturbation.  We extend the result to establish an oracle bound on the treatment effect.

We also offer several statistical contributions. We show that our constructed bases produce regression estimates that are comparable to cutting-edge machine learning methods.  We also implement a sparse Bayesian model with a prior constructed to achieve desirable predictive properties.  The  prior on the key tuning parameter grows  grows in the sample size and number of covariates so as to guarantee an oracle bound.  The prior also includes adaptive weights, similar to those in \citet{zouadaptive}, except they are estimated endogenously.  Furthermore, the prior is nonseparable, so a global tuning parameter helps the parameters adapted to the estimated level of overall sparsity in the model \citep[similar to the nonseparable priors in][]{bhattacharya2015dirichlet, Rockova17}.  In simulations we show that the gains over existing sparse Bayesian models \citep{cps2010horseshoe, polson2014bayesian} are sizable.


%The third major advance is theoretical.  Our key tuning parameter is set to grow at an oracle rate.  Estimating an outcome well is not precisely the same as predicting the impact of perturbing the treatment.  As one example, ensemble tree methods will fit the systematic component well, but since the fits are piecewise constant, they are not well-suited to estimating the effect of a local perturbation.  We extend the result to establish an oracle bound on the treatment effect.

%Statistical contributions of the paper benefit from tightly coupling the causal inference theoretical framework with estimation strategies estimation from machine learning. By combining particular classes of basis functions that are implied by causal inference theory and a recent sparse Bayesian model, we obtain oracle bounds on the key quantities of interest.  The model includes a prior such that the key tuning parameter grows  in the sample size and number of covariates so as to guarantee an oracle bound on perturbations of the treatment (i.e., quantity like a local derivative). \textbf{The theoretical connections are important because estimating an outcome well (as in much previous work in machine learning where oracle bounds are presented) is not precisely the same as predicting the impact of perturbing the treatment.}\footnote{As one example, ensemble tree methods will fit the systematic component well, but since the fits are piecewise constant, they are not well-suited to estimating the effect of a local perturbation.}  Finally, our sparse Bayesian model takes advantage of a range of several advances in the sparse modelling tradition. We use adaptive weights, similar to those in \citet{zouadaptive}, except that we estimate them endogenously to model.  Furthermore, the prior is nonseparable, so a global tuning parameter helps the parameters adapted to the estimated level of overall sparsity in the model \citep[similar to the nonseparable priors in ][]{bhattacharya2015dirichlet, Rockova17}. With these contributions in mind, we show that the proposed methods beats or is highly competitive with cutting-edge machine learning methods including the sparse Bayesian model class we draw from \citep{cps2010horseshoe, polson2014bayesian}.



We then extend the method to several structural models and treatment types.  In the case of a binary treatment, the proposed method outperforms existing cutting-edge machine learning methods on a benchmark dataset \citep{lalonde1986evaluating}.  In the instrumental variable setting, we show that the method returns  an observation level Wald estimate. Furthermore, we use the oracle bound constructively, rather than simply descriptively, in order to distinguish observations impacted by the instrument from those unaffected.  This distinction is crucial, as it is only off the former group that we can estimate a causal effect.    This oracle threshold serves as an improvement over existing high-dimensional instrumental variables methods, which do not distinguish between impacted and non-impacted observations.    We also extend the method of direct estimation to the mediation framework, where a treatment may have a direct effect on an outcome or through a post-treatment mediator.  We show that the method allows us to maintain a key identity in mediation analysis: that the mediation effect is the product of the treatment effect on mediator and mediator direct effect on outcome, a relationship otherwise lost in a nonlinear model.  Lastly, we show that the model extends naturally to the estimation of controlled direct effects via sequential g-estimation \citep{acharya2015explaining,vansteelandt2009estimating}.

Though most of the work focuses on estimating the local effect of a small perturbation of the treatment, we also extend the method to estimate the support of the treatment, giving the researcher a sense of over what range the model can be safely used to predict outcomes for a given observation, conditional on observed covariates.  We then use the estimated density of the treatment to estimate not the local effect of the treatment at the observed value but to the expected effect of the treatment over all treatment values, again conditional on covariates.


The paper proceeds as follows. Section~\ref{sec:framework}  introduces the proposed framework, estimation  strategy, and simulations.  Section~\ref{sec:extensions} shows how the method of direct estimation readily extends a binary treatment framework as well as structural models like instrumental variables, mediation, and sequential g-estimation. Section~\ref{sec:earlier} discusses the relationship between the proposed method and earlier work. We include several applications in Section~\ref{sec:app}. Section~\ref{sec:conclusion} concludes by showing how our approach connects disparate literatures and by discussing limitations and future extensions.

\begin{comment}

Recent work on drawing causal inference from observational and experimental data has emphasized two lessons.  First, causal effects are defined at the observation level.  Second, causal effects are defined in terms of ceteris paribus manipulations of a treatment.   As in most research designs we do not observe the outcome under counterfactual manipulations, a key observation-level quantity remains unobservable  \citep{Imbens2015, Rubin2005, holl:86}.

%The most commonly encountered estimands are aggregates of these individual effects.

Despite these insights, causal estimation has focused on estimating structural parameters that are  aggregates of these observation-level effects.  For example, inverse probability weights or subclassification can reduce confounding bias, but only an average effect for the sample is recovered \citep{rose:rubi:83, ho:imai:king:stua:07,  Iacus2011, Diamond2012, imai2012causal, Hirano2005}.  Estimating only an average effect ignores potentially informative underlying heterogeneity.  The shortcoming grows more pernicious in more complex designs.  In an instrumental variables analysis, for example, a causal effect is estimated but only for the unknown subset actually encouraged by the design \citep[e.g.][]{angrist1996identification}.   Uncertainty over the compliant subset, as well as unverifiable behavioral assumptions that no observation acts contrary to the encouragement, limits the method's internal and external validity \citep[e.g.][]{deat:10}.  Additional models for the treatment mechanism, inverse probability weights, and compliance status may help alleviate these problems \citep[e.g.][]{aronow2013beyond, Heckman07a, Heckman07b, Abadie2003, hira:imbe:rubi:zhou:00} but the models then introduce another layer of modeling assumptions and uncertainty.


Rather than target low-dimensional population-level parameters that correspond with an aggregate of individual effects, we propose a method whereby we use observed covariates to directly estimate observation-level quantities of interest. We call this the \emph{method of direct estimation}.  The estimates can be analyzed for heterogeneities or aggregated into coarser causal estimands, such as an average effect.    In more complex designs, the method offers deeper insight.  For example, in an encouragement design, the standard estimate is the Wald estimate, a ratio of two difference-in-means.  The sample differences in the numerator and denominator average over both those encouraged and those not encouraged.  Instead of a ratio of means, our estimate is a mean of observation-level ratios, allowing the researcher to estimate the subset of the sample actually encouraged.

The method extends directly to general treatment regimes.  In the case of a categorical or binary treatment, we estimate as the difference between the outcome under treatment and under some baseline level.  In the case of a continuous treatment, we estimate the partial derivative of the outcome with respect to the treatment, which is simply the analog of the sub-differential used with a the difference in the treatment levels taken to zero. We focus on the impact of a small ``perturbation'' as this measures the local marginal impact of the treatment on the outcome \citep[e.g.,][]{Hardle89} and  predicting near the observed data allays concerns about lack of common support and the dangers of extreme counterfactuals \citep{king2006dangers}.



A key goal of this project is to connect machine learning and causal inference at a theoretical level.  Rather than turn to a generic machine learning method to estimate causal estimands, we  use our identification assumptions to drive our modeling choices.  For example, we use the ignorability assumption from causal inference to motivate our choice of nonparametric basis.   We show that under this assumption, the partial derivative (or subderivative) of the outcome with respect to the treatment is only a function of pre-treatment covariates.  This suggests modeling the treatment in terms of a nonparametric basis that is locally linear, i.e. a degree 2 truncated-power basis spline or order 2 B-spline \citep[e.g.][ch. 5]{hast:tibs:frie:10}.   We are using the model to adjust for confounding, which is induced by a systematic interaction between the treatment and a prognostic variable.    We turn to a tensor-spline basis to account for these interactions.  As the covariates may also be interacting with each other, we fit a model with bases that capture $treatment \times covariate \times covariate$ interactions, as well as all lower-order terms.%, with ``interactions'' meant to be interpreted in the tensor-product functional sense.

As a second connection, we use the theory of LASSO estimation to derive a bound on the size of expected prediction error \citep{Buhlmann2013}.  We fit the model using a sparse Bayesian method \citep{RatkovicTingley17}, but we use this oracle inequality to drive our prior specification.  We show that doing so generates predictions notably better than existing sparse Bayesian methods  \citep{cps2010horseshoe, polson2014bayesian}.  We extend the theory, deriving an oracle bound not only for the prediction error but also for the  causal effect.  We show that solving a LASSO problem also bounds the error rate on the treatment effect.

Third, we then use this oracle bound constructively rather than simply descriptively.   We construct a feasible estimate for this bound, and then we show that the bound is a reliable means of differentiating compliers, non-compliers, and defiers in an instrumental variable analysis.  Doing so increases the internal validity of an instrumental variable analysis, as it helps identify off which observations a causal effect is identified.  It also aids  external validity, helping to establish what subgroups can be extrapolated to.

%Our conceptual and estimation framework offer several advances.  As we are simply predicting counterfactual outcomes, the method extends directly to general treatment regimes.  In the case of a categorical or binary treatment, we estimate as the difference between the outcome under treatment and under some baseline level.  In the case of a continuous treatment, we estimate the partial derivative of the outcome with respect to the treatment, which is simply the analog of the sub-differential used with a the difference in the treatment levels taken to zero. We focus on the impact of a small ``perturbation'' as this measures the local marginal impact of the treatment on the outcome \citep[e.g.,][]{Hardle89} and  predicting near the observed data allays concerns about lack of common support and the dangers of extreme counterfactuals \citep{king2006dangers}.




Our estimation strategy incorporates two innovations.  First is the construction of a nonparametric basis spanning a larger model space than  earlier work on multivariate spline models \citep{Stone1997,friedman1991multivariate}, as we include bases of multiple degrees and knot locations.  As the number of bases places us in an ``ultrahigh'' dimensional setting, we use a covariate screening method \citep{Fan:2008} to return a feasible subset of bases.  Second, as described above, we use a sparse Bayesian regression to estimate the model \citep{RatkovicTingley17}.  The regression model incorporates frequentist Oracle results to motivate hyperprior parameters.  The model also includes endogenously estimated adaptive weights \citep[as in][]{zouadaptive} and a  decay parameter that adapts these weights off the estimated global level of sparsity in the model \citep[similar to the nonseparable priors in ][]{bhattacharya2015dirichlet, Rockova17}.


%Our estimation strategy builds on earlier work on multivariate spline models \citep{Stone1997,friedman1991multivariate}, but incorporates several advances.  First, we consider a larger model space than earlier methods, allowing for several different sets of spline bases of different degree and knot location for each variable.  Since this places us in an ``ultrahigh'' dimensional setting, we use a covariate screening method \citep{Fan:2008} to return a feasible subset of bases.  The screen also allows us to consider the case where we have more covariates than observations. % Second, extend recent work \citep{RatkovicTingley17}, we implement a high-dimensional Bayesian model for estimating the basis weights.   Existing Bayesian methods set this key shrinkage hyperprior parameter as some constant \citep{cps2010horseshoe, Polson15, Rockova17} despite theoretical results that show how this parameter should grow in sample size and number of covariates \citep[e.g.][Sec 6.2]{Rockova17}.   We implement a prior constructed such that the theoretical rate is in fact achieved, and we show this leads to predictive gains relative to existing sparse Bayesian implemenations.   The prior is constructed such that estimates achieve the near-minimax ``oracle'' rate if the model is in-truth sparse.    If the model is in-truth dense, the weights adjust and the model still achieves the near-minimax rate.   A hyperprior parameter estimated from the data adjudicates between these two cases, generating a non-separable prior on the mean parameters, a technique found effective in other work \citep{bhattacharya2015dirichlet, Rockova17}.

%The model also estimates adaptive weights that, under certain conditions, achieve a faster convergence rate, a property we term ``super-persistent.''

Every method comes with shortcomings.  Estimating the outcome model without also modeling the treatment mechanism results in inefficient estimates \citep{robins1997toward}.  Yet, uncertainties in the treatment model may affect the estimates in nonlinear and unpredictable ways \citep{kang2007demystifying}.  We place a greater emphasis on model selection than efficiency.  We also note that our model requires assumptions on the differentiability' of the conditional mean function.  In our framework, these assumptions are unavoidable, but we work with a dense and varied set of discontinuous bases to accommodate a wide class of functional forms.  Our simulations include data generated from models both in and out of our model space, and we find the proposed method works well.

We illustrate the method through simulation studies and two applied examples.   We illustrate in a simulation study how our proposed method generally outperforms a variety of cutting-edge machine learning and high-dimensional methods in predicting both fitted values and the marginal responsiveness of the outcome to the treatment.  We then illustrate the method using a binary treatment and then continuous treatment.  For the binary treatment, we turn to a benchmark dataset in the matching literature but comes from important work in economics \citep{lalonde1986evaluating}.  We show that the method of direct estimation outperforms existing measures in recovering an average treatment effect from an observational dataset.  As we focus on individual-level estimates, we also show that our method performs well in recovering individual-level estimates for a held-out subsample of the data. We next turn to applying the method to an instrumental variables analysis. We show how the method relates to earlier approaches to instrumental variables that focus on structural parameters, evaluate performance of the method against two-stage least squares, and apply the method to recent work on the relationship between education and political participation \citep{larreguy2017}.  We show how the proposed method uncovers individual compliance status as well as an understanding of how causal effects vary with the levels of the instrument. Finally, we analytically show how the method extends to mediation analysis \citep{baron1986moderator,vanderweele2015explanation}.




The paper proceeds as follows. Section~\ref{sec:framework}  introduces the proposed framework, estimation  strategy, and simulations.  Section~\ref{sec:extensions} shows how the method of direct estimation readily extends a binary treatment framework as well as structural models like instrumental variables, mediation, and sequential g-estimation. We include several applications. Section~\ref{sec:earlier} discusses the relationship between the proposed method and earlier work. Section~\ref{sec:conclusion} concludes by showing how our approach helps to connects disparate literatures and by discussing limitations and future extensions.

\end{comment}





\section{The Proposed Method}
\label{sec:framework}

We first describe the setup for the proposed method. We next describe each step of the estimation strategy: generating bases, screening, our sparse regression model, and predicting counterfactuals.  We then present our oracle bounds on both the fitted values and the treatment effect.



\begin{comment}

but by looking at the observation-level ratio we are able to directly account for heterogeneity in the instruments impact on the treatment.


Our argument follows in two parts.

Given that we target this manipulation, we introduce a nonparametric regression method for estimating

With this simple case in hand, we then


Marc please focus on the below{\bf Dustin--think we did this last night, right?}
%CATE versus local perturbation

Local perturbations and lots of interactions

Local perturbation <---> CATE

Local do operation <-->  pointwise CATE <--> local perturbation  <--> partial derivative of conditional mean

Doing CATE--super local with perturbations

Local CATE

Flexibility--local cate estimates local causal effect--we get closer to conditional mean than existing methods


Saturated parameter space crossed with local treatment perturbation

Why we need crazy bases

At each point, we look locally and it linearizes the problem

\end{comment}



\subsection{Setup}

Assume a simple random sample where observation $i\in \{1, 2, \ldots, n\}$ possesses a potential outcome function, $Y_i(T_i)$, which maps treatment level $T_i$ to outcome $Y_i(T_i)$.   The treatment $T_i$ is a realization  of $\widetilde T_i \sim F_{X_i}$, such that $\widetilde T_i$ has continuous support and its distribution may vary with  $X_i$, a vector of $p$ observation-level fixed covariates.    There may be more covariates than observations, so $p$ may be greater than $n$, and the covariates may include moderators, risk factors, prognostic variables, or simply superfluous variables; the important point is that they are fixed and causally prior to $\widetilde T_i$.  The observed data are $\l\{Y_i, T_i, X_i^\top\r\}^\top$ where $Y_i = Y_i\l( T_i\r)$


We will denote the causal effect of moving from one value of the treatment, $T_i^\prime$, to another, $T_i^{\prime\prime}$,  as
\begin{align}
\Delta_i(T_i^{\prime\prime}, T_i^\prime)=Y_i(T_i^{\prime\prime}) - Y_i(T_i^\prime).
\end{align}

To identify our causal effects, we make three standard assumptions.  First, we assume the potential outcomes are ignorable given $X_i$: $Y_i(\widetilde T_i)\indep \widetilde T_i|X_i, \widetilde T_i \in\{T_i^{\prime\prime},T_i^\prime\}$.  Second, we assume the counterfactual is possible, $\Pr(\widetilde T_i=T_i|X_i, T_i \in\{T_i^{\prime\prime},T_i^\prime\})>0$.  Third, we make stable unit and treatment value assumptions: potential outcomes do not depend on treatment or outcomes of others and there is a single version of each treatment level.


The quantity $\Delta_i(T_i^{\prime\prime}, T_i^\prime)$ poses two challenges.  First, as we can observe at most one potential outcome per observation,  at least one of the potential outcomes remain hidden.  This requires us to instead condition on the covariates, estimating the conditional effect for observation $i$ as $\E(\Delta_i(T_i^{\prime\prime}, T_i^\prime)|X_i)$, the expected change in outcome for an observation with profile $X_i$.  A second concern is whether $T_i^\prime$ and $T_i^{\prime\prime}$ might occur with positive probability.  To reduce concerns over common support,  we first focus on the partial derivative $\lim_{\delta \rightarrow 0}\nabla_{T_i}(\delta_i)$, where
\begin{align}
\nabla_{T_i}(\delta_i) &= \frac{1}{\delta_i} \E\l\{ Y_i\l( T_i +\delta_i\r) - Y_i\l( T_i \r) |X_i\r\},
\end{align}
which is a difference off the observed $T_i$, which must have positive support. The estimand $\lim_{\delta \rightarrow 0}\nabla_{T_i}(\delta_i)$ is the local effect of a  perturbation of $T_i$ on $X_i$.  For example, researchers are used to fitting a model of the form
\begin{align}
Y_i = \mu_i + \beta T_i +X_i^\top \gamma + \epsilon_i
\end{align}
where $\beta$, a global slope term, is interpreted as the marginal effect of $T_i$ on $Y_i$.  We conceive of the outcome in the following form:
\begin{align}
Y_i = \mu_i + \beta_i T_i + g(X_i) + \epsilon_i; \label{eq:outcome1}\\
\beta_i =f(X_i)= \lim_{\delta \rightarrow 0}\nabla_{T_i}(\delta_i) \label{eq:outcome2}
\end{align}
and $f(X_i), g(X_i)$  nonparametric functions of the covariates.  That the local effect is only a function of $X_i$ follows directly from our ignorability assumption:
\begin{align}
&Y_i(\widetilde T_i)\indep \widetilde T_i|X_i, \widetilde T_i \in\{T_i+\delta_i,T_i\}\\
&\Rightarrow \lim_{\delta_i\rightarrow 0}\E\l\{ \frac{Y_i(T_i+\delta_i)-Y_i(T_i)}{\delta_i}\bigg|X_i\r\}=\\
&\quad\quad\lim_{\delta_i \rightarrow 0}\nabla_{T_i}(\delta_i)=f(X_i)
\end{align}
where $f(X_i)$ is not a function of $T_i$.  Taking an antiderivative gives the model in Equations \ref{eq:outcome1}-\ref{eq:outcome2}.

We focus on this local effect through the first part of the paper.  Doing so helps us forestall issues related to extrapolating off the support of $\widetilde T_i,$ as we are staying infinitesimally close to the observed $T_i$. While we return to the issue later, we first focus on developing a rich regression model connecting the treatment to the outcome, locally at the observed value.   Confounding bias is introduced when the treatment and covariates interact but the interaction is unmodeled. Failure to model this interaction is even more pernicious when the interaction is non-linear.  In order to more reliably estimate the treatment effect, we allow for the local interaction between the treatment and a rich set of covariate spline bases.  Addressing this local confounding returns a local treatment effect estimate--a conditional effect for that observation, or, equivalently, the impact of a tiny ``do'' operation or the partial derivative of the response function with respect to the treatment.\footnote{In Pearl's ``do'' notation, $\nabla_{T_i}(\delta_i)  =\frac{1}{\delta_i} \l\{ \E(Y_i|do(T_i+\delta_i), X_i)-\E(Y_i|do(T_i), X_i)\r\}$.}  Our emphasis is on implementing a sufficiently rich regression model to estimate and remove confounding bias at each observation.

Once we have a means of estimating this local effect, we illustrate several useful ways in which we can extend it to other causal settings.  In an instrumental variable setting, we can model first-stage heterogeneity in the encouragement, allowing us to relax the no defiers assumption.  In a mediation setting, we show that we can maintain a key identity--that the mediation effect is the product of the treatment effect on mediator and mediator direct effect on outcome--that is otherwise lost in nonlinear models. We also show that the method can estimate controlled direct effects in a single step, but can also extend to a broad class of structural mean models through the use of sequential $g$-estimation.

We then move past working with partial derivatives to estimating and averaging over the support of $\widetilde T_i$.  We use a bootstrap-estimated interval to characterize a range over which the treatment effect can be estimated.  We also use the bootstrap to move from estimating the effect of the treatment at the observed value, $\lim_{\delta_i \rightarrow 0}\nabla_{T_i}(\delta_i)$ to the expected effect of $T_i$, $\E_{\widetilde T_i}\l(\lim_{\delta_i \rightarrow 0}\nabla_{T_i}(\delta_i)\r)$. Before exploring these extensions we turn to our estimation procedure.


%Once we have this basic machinery in place, we generalize in several directions in Section XX. First, we show that the estimand $\nabla_{T_i}(\delta_i)$ can be extended to instrumental variable models and mediation models, providing insight to the workings of each.  Second, we extend the method to categorical treatments.  Third, we show how to estimate the support of $\widetilde T_i$, giving a range over which extrapolation is reasonable.  Lastly, we give a means of moving from estimating the partial derivative at the observed value, $\lim_{\delta \rightarrow 0}\nabla_{T_i}(\delta_i)$ to the expected effect of $T_i$, $\E_{\widetilde T_i}\l(\lim_{\delta \rightarrow 0}\nabla_{T_i}(\delta_i)\r)$.


\subsection{Estimation}
\label{sec:estimation}

We fit a model using a high-dimensional, tensor-spline nonparametric regression.  Estimation proceeds in four steps.  First, we generate a large (on the order of millions) set of nonparametric bases.  Second, we use a marginal correlation screen to select a subset of bases, reducing their number from millions to hundreds or thousands.  Third, we use a sparse Bayesian regression to fit the model.  Fourth, we use the estimated model to predict the counterfactual quantities of interest.

%The advances in this approach are twofold.  First, we are the first to use a variable screen in a high-dimensional, tensor-spline regression.  The spline basis, even when fit using a vanilla cross-validated LASSO,  outperforms many ``black-box,'' ensemble methods in our simulations and applications.  Second, we implement a sparse Bayesian prior calibrated to achieve an Oracle bound.  While our sparse Bayesian prior is similar to existing variants, we are the first to construct a prior such that the key tuning parameters ensure desirable Oracle prediction properties.

%We now present our statistical model for the observed outcome and then move through our estimation procedure.
%\subsubsection{The Model}
We assume the data are generated as
\begin{align}
Y_i = \mu_Y +  R^o(T_i, X_i)^\top c^o +  \epsilon_i \label{eq:out1}
\end{align}
with $\epsilon_i \iid \mathcal N(0, \sigma^2)$, where $R^o(T_i,X_i)$ is a possibly infinite set of basis functions.  We divide these into three groups:  those that survive our covariate screen (described below),  $R(T_i, X_i)$; those that we consider in the screen, $R_\infty(T_i, X_i)$; and those that are in the data-generation process but not in our considered model space, $R_\perp(T_i, X_i)$.  Each will later be crucial components of our bound analysis.

\subsubsection{Generating Bases}
In constructing the set of bases that we will use to model the conditional mean, $R_\infty(T_i, X_i)$, we turn to a tensor-product of spline bases. Denote $\{1,b(X_{ik}, \kappa_k)\}_{\kappa_k \in \mathcal K_k}$ as an intercept term and a set of $|\mathcal K_k|$ mean-zero nonparametric bases for variable $X_k$ evaluated for observation $i$ at the knots in $\mathcal K_k$.  Similarly, for the treatment, denote $\{1, b_T(T_{i}, \kappa_T) \}_{\kappa_T \in \mathcal K_T}$ as an intercept and $|\mathcal K_T|$ bases evaluated for the treatment observed for observation $i$ evaluated at knots in $\mathcal K_T$.  We construct the mean vector as
\begin{align}
R_\infty(T_i, X_i) = \{1, b_T(T_{i}, \kappa_T) \}_{\kappa_T \in \mathcal K_T} \otimes \{1, b(X_{ik}, \kappa_k)\}_{\kappa_k \in \mathcal K_k}\otimes\{1, b(X_{ik^\prime}, \kappa_{k^\prime})\}_{\kappa_{k^\prime} \in \mathcal K_{k^\prime}}\label{eq:tensor1}
\end{align}
where $\otimes$ denotes the tensor product.  In effect, we construct a basis for the mean that contains nonparametric bases in the treatment and covariates while allowing for $ treatment \times covariate$, $covariate \times covariate$, and $treatment \times covariate \times covariate$ interactions.


Our local-slope formulation implies then that the outcome is locally linear at each $T_i$ with local effect given by an arbitrary function of $X_i$.  We therefore model the treatment using a degree 2 thresholded power-basis  and a degree 2 $B$-spline.  The first is a set of hinge functions radiating off knots, the second is an upside-down $V$ between two knots, and zero elsewhere.  Both sets of bases satisfy the local linearity, unlike bases with more than one nonzero derivative in $T_i$ (e.g. gaussian radial basis functions, higher degree $B$-splines).

For each covariate, we use a $B$-spline basis of varying knots and degree \citep{Eilers96, deBoor78}.  The $B$-spline basis offers a bounded basis for nonparametric modeling, with well-studied optimality properties \citep{Gyorfi2002}.  The $B$-spline basis requires selecting both degree and knot locations, and as we do not know the best choice, we err in favor of selecting too many basis functions.  Specifically, for $k\in\{3,5,7,9\}$, we model each confounder using $k$ centered B-spline bases of degree $k$ with knots at every $100/(1+k)^{th}$  percentile of the covariate, generating $24=3+5+7+9$ bases.


\subsubsection{Screening Bases}

We  combine the treatment and covariate bases using the tensor product as given above.    At our default, we take $28$ bases for  the treatment and $25$ for each covariate.  This gives us $(1+28)\times(1+24\times p)\times(1  + 24 \times p)$ total basis.  Even a modest $p$ generates a large number of bases; $p=10$ gives over $1,680,000$ bases and $p=20$ over $6,700,000$.  This places us in an ``ultrahigh'' dimensional setting so we implement a Sure Independence Screen (SIS) \citep{Fan:2008}.\begin{footnote}{We implement a SIS rather than the grouped SIS strategies of \citep{fan2012nonparametric,fan2014nonparametric}, as we are not willing to assume the group structure necessary to bring in entire sets of bases, and the grouped SIS has not been extended to tensor product spaces.}\end{footnote}   The SIS occurs in two steps.  First, we construct all of the tensor product bases in Equation \ref{eq:tensor1}. Second, we sort the bases by their absolute correlation with the outcome.  We want to select a sufficiently large and rich dictionary of nonparametric bases to approximate a wide set of models; in practice and all of the analyses below, we select $100\times (1+n^{1/5})$ bases that have the largest unconditional correlation with the outcome.\begin{footnote}{Memory and computational speed are a concern.  To preserve memory, we  construct millions of bases, but at any given point in their construction we only save the $100\times (1+n^{1/5})$ bases with the largest absolute correlation with the outcome.  Construction of the tensor basis is done in C++ via {\tt Rcpp}.   }\end{footnote}

\subsubsection{A Sparse Bayesian Model}
Even after screening, we are left with hundreds of possible nonparametric bases.  In order to enforce some form of regularization, we implement a sparse Bayesian prior.  We estimate our treatment effects by adapting  LASSOplus, a sparse Bayesian prior we implemented in earlier work \citep{RatkovicTingley17}, to the nonparametric prediction problem.


The likelihood component is a normal regression model.  The prior structure has three properties. First, it ensures that estimates satisfy an oracle bound on the prediction error.  Second, a set of inverse weights on the mean parameters regularize large effects less aggressively than small effects.  Third, a global sparsity parameter produces a non-separable prior, allowing for the estimates to adapt to the overall level of shrinkage in the model.

The hierarchy is given as
\begin{align}
Y_i |R(T_i, X_i), c, \sigma^2 &\sim \mathcal N(R(T_i, X_i)^\top c, \sigma^2)\\
c_k| \lambda, w_k, \sigma &\sim DE\l(\lambda w_k/\sigma\r)\\
\lambda^2|n, p, \rho& \sim \Gamma\l( n\times(\{log(n)+2\log(p))-p, \rho\r) \\
w_k|\gamma &\sim \textrm{generalizedGamma}(1,1,\gamma)\\
\gamma &\sim \exp(1) \label{eq:gammaprior}
\end{align}
where $DE(a)$ denotes the double exponential density, $\Gamma(a, b)$ denotes the Gamma distribution with shape $a$ and rate $b$, and $\textrm{generalizedGamma}(a,b,c)$ the generalized Gamma density $f(x;a,d,p )=\frac{p/a^d}{\Gamma(d/p)}x^{d-1}\exp\{-(x/a)^p\}$.  We have two remaining prior parameters.  We take a Jeffrey's prior $\Pr(\sigma^2) \propto 1/\sigma^2$ on the error variance and we set $\rho=1$ in the generalized Gamma density.   We fit the model using an EM algorithm and use hats to denote the fitted values, i.e. $\widehat c$ is the EM estimate for $c$.

We summarize several properties of our model and present formal derivations in  Appendix \ref{app:lassoplus}. First, the prior on $\widehat \lambda^2$ is scaled in $n,p$ such that the estimate $\widehat \lambda$ achieves the oracle growth rate of $\sqrt{n \log(p)}$ {\it ex post} when $p$ is of order $n^\alpha$, $\alpha>0$, as in the nonparametric setting here.  Second, the rate at which the oracle bound holds is controlled by $1-exp\{-\sqrt{\log(n)}\}$, which approaches $1$ in the limit.\begin{footnote}{The  prior is proper whenever $n\times(\log(n)+2\log(p))-p>0$.  For example, with $n\in\{100; 1,000;10,000\}$, this requires $p<1,978; 27,339; 347,529$ respectively, which is sower than the LASSO rate of $n \sim \log(p)$, but clearly allows for $p>n$.  This slower growth rate does suggest the utility of the an initial screen for covariates.  Last, even with an improper prior, the posterior will be proper, so estimation can still proceed.}\end{footnote}

The prior weights $w_k$ are constructed so that the MAP estimate is similar to the  adaptive LASSO of  \citet{zouadaptive}.   Each mean parameter, $c_k$, will have its own adaptive penalty $\frac{\widehat \lambda\widehat w_k}{n\widehat \sigma}$.The estimated weights $\widehat w_k$ are inversely related to the magnitude of the parameter estimates $\widehat c_k$.  We prove that the adaptive penalty term is of order  $\sqrt{ \log(n)/n}$ when  $\widehat c_k \rightarrow 0$, which is not asymptotically negligible.  On the other hand, the adaptive penalty is of order $1/n$ when $\widehat c_k \not \rightarrow 0$, and hence is asymptotically negligible.

Third, the parameter $\widehat \gamma$ adapts to the global sparsity level of the data.  If we take $\gamma \rightarrow 0$, then the prior approaches a degenerate `spike-and-slab' prior uniform over the real number line but an infinite point mass at 0.  In this scenario, we are not shrinking any effects except for those exactly zero.   At the other extreme, $\gamma \rightarrow \infty$, the prior approaches a Bayesian LASSO of \citep{park2008bayesian}, regularizing every term. The global tuning parameter, $\widehat \gamma$, is estimated from the data and adjudicates between these two extremes; the utility of nonseparable priors over the tuning parameters have also been discussed by \citet{bhattacharya2015dirichlet, Rockova17}, though we note that these priors were not tuned to achieve the oracle property (or similar concentration property) {\it ex post}.  Details are provided in Appendix \ref{app:lassoplus}.


\subsubsection{Counterfactual Prediction}

Given fitted values $\widehat c$, we can write
\begin{align}
\widehat Y_i = \widehat \mu_Y+R(T_i, X_i)^\top \widehat c
\end{align}
with the intercept chosen as $\widehat \mu_Y =\frac{1}{n} \sum_{i=1}^n  \l(Y_i - R(T_i, X_i)^\top \widehat c\r)$.  We then estimate
\begin{align}
\widehat \nabla_{T_i}(\delta_i) =\frac{1}{\delta_i} \l\{\big(R(T_i+\delta_i, X_i)- R(T_i, X_i)\big)^\top \widehat c\r\}
\end{align}
and we approximate the derivative by choosing $\delta_i$ close to zero; we take $\delta_i = 10^{-5}$ in practice.  The sample average marginal causal effect can be found through averaging over the sample, $\frac{1}{n} \sum_{i=1}^n \widehat \nabla_{T_i}(\delta_i)$.

\subsection{Theoretical Results}

We provide two bounds on the risk of our estimator.  First, we provide a bound on the predicted values.  Second, we provide a bound on the effect estimate.  The first follows almost directly from our prior structure, which guarantees that the key tuning parameter grows at the Oracle rate. The second bound shows that the model  will do well in terms of predicting the local treatment effect as well.



\paragraph{An Oracle Inequality}   Denote as $R(T_\cdot, X_\cdot)$ the $n\times p$ matrix of bases and $\widehat \delta=\widehat c- c$.    Since our tuning parameter grows at the oracle rate $\sqrt{n \log(p)}$, we can bound the excess risk as
\begin{align}
%\frac{1}{N}\l\{||R(T_i, X_i) \widehat \delta||_2^2 + \lambda ||\widehat W \widehat \delta||_1^1\r\} \le \frac{C \widehat \sigma^2 \overline \gamma^2 |S|}{N^2 \phi_0^2}\label{eq:oracle1}
\begin{split}
&\frac{1}{n}\l\{||R(T_\cdot,X_\cdot) \widehat \delta||_2^2+ \lambda||\sum_{k=1}^p \widehat w_k \widehat \delta_k||_1^1\r\}\le\\
&\quad C\frac{\widehat \lambda^2 \widehat \sigma^2 \overline \gamma^2 |S|}{n^2 \phi_0^2}+ C_\infty \frac{||R^o_{\infty/n}(T_\cdot,X_\cdot)c_{\infty/n}||^2_\infty}{n} +C_\perp\frac{  ||R^o_{\perp}(T_\cdot,X_\cdot)c_{\perp}||_\infty^2 }{n}  \label{eq:oracle1}
\end{split}
\end{align}
with a probability at least $\exp\{-\sqrt{\log(n)}\}$.  The bound splits into three components, a bound off the post-screened basis $R(T_\cdot,X_\cdot)$, a bound attributable to bases that did not survive the screen but would as $n$ grows, $R^o_{\infty/n}(T_\cdot,X_\cdot)$, and a bound attributable to the portion of the model that falls outside the span of the basis used in the screen, $R^o_{\perp}(T_\cdot,X_\cdot)$.  In the first component,
 $\phi_0$ is the smallest eigenvalue of the Gram matrix of the submatrix of $R(T_i,X_i) \widehat W^{-1/2}$, $|S|$ is the number of non-zero elements of $c$ and we use $C_\cdot$ to denote an unknown constant that does not grow in $n,p,S$.  For details, see \citet{RatkovicTingley17}.  The next two components are attributable to differences between $R^o(T_\cdot, X_\cdot)$ and $R(T_\cdot,X_\cdot)$.   The second term is the same order of the first--both are of order  $O(\log(n)/n)$; see Appendix \ref{app:proofs} and \citet{Fan:2008}.  The third term is irreducible and of order $O(1)$, corresponding with inescapable misspecification error.    We minimize our concerns over the third term through selecting a sufficiently dense set of bases.

\paragraph{An Oracle Inequality on  $\nabla_{T_i}(\delta_i)$}  We are interested in not only bounding the prediction risk, as above, but also the error on the $\nabla_{T_i}(\delta_i)$.      To do so, we decompose our bases into two components, one submatrix such that no element covaries with the treatment and one submatrix where for each basis at least one element covaries with the treatment:
\begin{align}
R(T_\cdot, X_\cdot) = [R^X(X_\cdot): R^{TX}(T_\cdot, X_\cdot)].
\end{align}
 Denote $\widehat Y_i = [R^X(X_i): R^{TX}(T_i, X_i)]^\top \widehat c$, and $\widehat c^X$ and $\widehat c^{TX}$ the subvectors of $\widehat c$ associated with $R^X(X_i)$ and $ R^{TX}(T_i, X_i)$.  We denote as $\Delta R(T_\cdot, X_\cdot)=[\Delta R_{ij}]=[\partial R_{ij}/\partial T_i]$ the elementwise partial derivatve of $R(T_\cdot, X_\cdot)$ with respect to the treatment and note $\widehat \nabla_{T_i}(\delta_i) = \Delta R_i(T_i, X_i)^\top \widehat c^{TX}=\Delta R^{TX}_i(T_i, X_i)^\top \widehat c^{TX}$.   Denote as $\widehat c^{\widehat S,TX}$  the subvector of $\widehat c^{TX}$ selected in the outcome model.


We show in Appendix \ref{app:proofs} that the mean parameter estimates $\widehat c^{\widehat S, TX}$ are actually a solution to the following LASSO problem:
\begin{align}
\widehat c^{S,TX} = \argmin_{c^{\widehat S,TX}} \big|\big|\widetilde{ \Delta Y} - \Delta R^{\widehat S, TX} c^{\widehat S,TX}  \big|\big|_2^2 + \big|\big| \widetilde W^{\widehat S, TX}c^{\widehat S,TX} \big| \big| \label{eq:lassoderiv}
\end{align}
where $\widetilde{\Delta Y}$ is a pseudo-response generated from the fitted derivative and estimated residuals from the outcome model, $\widehat \epsilon_i$, as
\begin{align}
\widetilde{\Delta Y_i} =\Delta R_i(T_i,X_i)^\top \widehat c + \widehat \epsilon_i.
\end{align}
The weight matrix $\widetilde W^{\widehat S, TX}$ is diagonal with entries $\widetilde W_{k}^{\widehat S, TX} = \big|\sum_{i=1}^N R_{ik}^{\widehat S, TX}(T_i, X_i)  \frac{\partial \widehat \epsilon_i}{
\partial T_i} \big|$.  For intuition, note that $\partial \widehat \epsilon_i/\partial T_i$ is the causal effect of the treatment on the prediction error, $R_i(T_i, X_i)(c-\widehat c)$.  The less correlated the basis $R_k^{\widehat S, TX}(T_i, X_i)$ with the impact of a treatment perturbation on prediction error, the less that basis's coefficient is penalized in predicting the treatment effect.  The more correlated the basis with the sensitivity of prediction error on perturbing the treatment, the more penalized the coefficient on that basis.

As the coefficients are simultaneously minimizing a LASSO problem on the first derivative, we can
establish the oracle bound
\begin{align}
\begin{split}
&\frac{1}{n}\l\{||\Delta R^{\widehat S, TX} (T_\cdot, X_\cdot)(\widehat c^{\widehat S, TX}-c^{\widehat S, TX})||_2^2 +  ||\widetilde W^{\widehat S, TX} (\widehat c^{\widehat S, TX}_k-c^{\widehat S, TX}_k)||_1^1\r\} \le \\
&C_\Delta  \frac{\widehat \sigma^2 \overline \gamma_\Delta^2 |S_\Delta|}{N^2 \phi_\Delta^2}+ C_\infty \frac{||\Delta R^o_{\infty/\widehat S}(T_i,X_i)c_{\infty/n}||^2_\infty}{n} +C_\perp\frac{  ||\Delta R^o_{\perp}(T_i,X_i)c_{\perp}||_\infty^2 }{n}
\label{eq:oracle2}
\end{split}
\end{align}
where all the constants in the inequality are analagous to those in Inequality \ref{eq:oracle2} but defined in terms of the design matrix that parameterize the derivative.

We last note that the model is similar to the relaxed LASSO of \citet{Meinshausen07}, differing in that rather than re-fitting the outcome to first-stage LASSO selected variables, we fit them to the pseudo-observation above.  In this formulation, we are conditioning on the fitted residuals, which suggests the utility of bagging to improve estimation--a point we return to below.

This oracle bound gives us a bound not on the predicted values but on the first derivative.  Bounds of this type have been used for primarily descriptive purposes, to establish rates of convergence of different estimates.  We use the bound constructively, to help us differentiate values for which the local treatment effect is zero from those where it is not.  The bound helps us in the instrumental variable setting to estimate off which observations we are identifying an effect.  We introduce a feasible estimate below, after we have established the reliability of the proposed method for point estimation through a simulation study.



\subsection{Simulations}
\label{sec:simulation}

%\subsubsection{Setup}

We next present simulation evidence illustrating the MDE's utility in estimating treatment effect.  We include four sets of simulations presented in increasing complexity, a linear model, a low-dimensional interactive model, a nonlinear model, and a model with interactions and discontinuous breaks, respectively:
\begin{align}
\textrm{Linear: }&\;Y_i=T_i+\sum_{k=5}^8 X_{ik} \beta_k+\epsilon_i\\
\textrm{Interactive: }&\;Y_i=T_i-T_i \times X_{i3}+\sum_{k=5}^8 X_{ik} \beta_k+X_{i1}X_{i2}+\epsilon_i\\
\begin{split}
\textrm{Nonlinear: }&\; Y_i=20\sin\l(\l(X_{i1}-\frac{1}{2}\r)\frac{T_i}{20}\r)+\frac{1}{2}\times (1+|X_{i3}\times  X_{i4}|_+)\times(1+ T_i)+\\
&\quad\quad\quad\quad\sum_{k=5}^8 X_{ik}\beta_k+\epsilon_i
\end{split}\\
\textrm{Discontinuous}&\;Y_i=\l(1+|T_i |\r)\times \mathbf 1(|T_i|>1/2)\times(X_{i5}+1) +\sum_{k=5}^8 X_{ik} \beta_k+X_{i1}X_{i2}+\epsilon_i
%        (1+abs(treat)*(abs(treat)>.5))*(X[,5]+1)+X[,1:4]%*%ran.beta[1:4]+X[,3]*X[,4]
\end{align}
where the the error is independent, identical Gaussian such that the true $R^2$ in the outcome model is $0.5$.  All covariates are from a multivariate normal with variance one and covariance of $0.5$ between all pairs and the elements of $\beta$ are drawn as independent standard normal.  The treatment is generated as
\begin{align}
T_i =-3+(X_{i1}+X_{i4})^2+\epsilon_i^T; \;\;\epsilon_i^T\iid \mathcal N\l(0, 4\r)
\end{align}
We consider $n \in \{100, 250, 500, 1000, 2000\}$ and $p \in \{10, 25, 50, 100\}$.  In each case, the true data is generated from the treatment and the first 8 covariates, but the dimensionality of the underlying model may be much higher.  For example, the smooth curve in simulation 4 is infinite dimensional (a sine curve).


We asses methods across two dimensions, each commensurate with our two estimation contributions: the nonparametric tensor basis for causal estimation and the sparse Bayesian method.  We want to assess, first, how well sparse regression methods perform given the same covariate basis.  We include in this assessment the  cross-validated LASSO, as a baseline, as well as two sparse Bayesian priors: the horseshoe and Bayesian Bridge \citep{cps2010horseshoe, polson2014bayesian}, as well as LASSOplus, described above.  Each of these methods are handed the same nonparametric basis created in our pre-processing step.  We compare these methods to regression-methods that generate their own bases internally,   kernel regularized least squares \citep[KRLS, ][]{Hainmueller13}, POLYMARS \citep{Stone1997}, and Sparse Additive Models \citep[SAM, ][]{Ravikumar09}.  These methods are simply given the treatment and covariates, not our nonparametric basis.  The last comparison set are non-regression methods, bayesian additive regression trees (BART, \citep{chipman2010bart}), gradient boosted trees \citep[GBM, ][]{Ridgeway99}, and the SuperLearner  \citep{polleysuperlearner}.  For the SuperLearner we include as constituent methods random forests, the LASSO, POLYMARS, and a generalized linear model.


\subsubsection{Results}

We report results on each methods' ability to recover the conditional mean $\mu_i =\E(Y_i|X_i, T_i)$ and partial derivative  $\partial \mu_i/\partial T_i$.  To summarize the performance, we calculate three statistics measuring error and bias.  The statistics are constructed such that they are $0$ when $\mu_i=\widehat \mu_i$ and $\partial \mu_i/\partial T_i=\widehat \partial \mu_i/\partial T_i$ and take a value of $1$ when $\widehat \mu_i = \frac{1}{n} \sum_{i=1}^n Y_i = \overline Y_i$.

\begin{align}
\textrm{RMSE}:\;&\; \sqrt{\frac{ \sum_{i=1}^n (\mu_i - \widehat Y_i)^2  }{\sum_{i=1}^n (\mu_i- \overline Y_i)^2 }}\\
\textrm{RMSE}_\nabla:\;&\; \sqrt{\frac{\sum_{i=1}^n \l(\frac{\partial \mu_i}{\partial T_i} - \widehat{\frac{\partial Y_i}{\partial T_i}}\r)^2 }{\sum_{i=1}^n  \l\{\frac{\partial \mu_i}{\partial T_i}\r\}^2 }}\\
\textrm{Bias}_\nabla:\;&\;\frac{ \bigg|\sum_{i=1}^n \frac{\partial \mu_i}{\partial T_i} - \widehat{\frac{\partial \mu_i}{\partial T_i}}  \bigg|}{\sqrt{\sum_{i=1}^n  \l\{\frac{\partial \mu_i}{\partial T_i}\r\}^2 }}
\end{align}



\begin{figure}[t!]
%\vspace{-.9in}
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{RMSE_Fitted.pdf}
    \vspace{-.25in}
  \end{center}\spacingset{1}
  \caption{{\bf A Comparison of RMSE for Fitted Values.}      } \label{fig:fd_fitted}
\end{figure}

In the figures, a value of $1$ can be interpreted as performing worse than the sample mean, in either a mean-square or bias sense, and values closer to 0 as being closer to the truth.  A value above $1$ when estimating the derivative means the method did worse than simply estimating $0$ for each value, the first derivative of the null model $\widehat \mu_i=\overline Y_i$.  For presentational clarity we remove any results greater than 1.5 for the RMSE results and .2 for the bias results.

Figure~\ref{fig:fd_fitted} reports the RMSE for the fitted values, by method.  We can see that MDE, the proposed method, performs well in each situation.  The first column contains results from the simple additive linear model, where MDE, POLYMARS, and the SuperLearner are all competitive.  We suspect SuperLearner is competitive here because the true model is in the model space for OLS, one of the components included in the SuperLearner.   MDE, POLYMARS and LASSO are nearly tied for the best performance in the interactive and non-linear settings. Finally while MDE is in a top performing set for the discontinuous case, POLYMARS had a slight advantage.



Figure~\ref{fig:fd_all} presents results on each methods' error in estimating the derivative.  In the linear model, we are consistently dominated by POLYMARS and beat SuperLearners as $p$ grows.   In the interactive and non-linear settings, MDE is again aligned with the LASSO and POLYMARS, while POLYMARS performs worse than the null model in the non-linear setting, with the horseshoe competitive.  POLYMARS performs the best in the discontinuous case but MDE was not far behind.  We find that, across settings, LASSOplus is the only method that generally performs first or second best in RMSE, though POLYMARS, particularly in the discontinous case, and cross-validated LASSO are reasonable alternatives.

\begin{figure}[t!]
%\vspace{-.9in}
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{RMSE_FirstDerivative_all.pdf}
    \vspace{-.25in}
  \end{center}\spacingset{1}
  \caption{{\bf A Comparison of RMSE for Unit Level Derivatives Across Methods.}      } \label{fig:fd_all}
\end{figure}


This highlights how  estimating fitted values well need not result in recovering the treatment effect well.  Most of the systematic variance in our simulations is attributable to the background covariates, and a method could perform well simply by getting these effects correct but missing the impact of the treatment.  For example, BART and Sparse Additive Models perform relatively well in estimating fitted values in the interactive, and nonlinear but not much better than the null model in estimating the partial derivative.

We perform well in our simulations, and we share the concern that this performance is the result of favorable decisions we made in the simulation design. To help alleviate these concerns, we use ordinary least squares to decompose the first derivative into two different components, one that correlates with the treatment and one that does not. Examining the two separate components helps identify the extent to which each methods captures the low-dimensional linear trend in the first dimension and the extent to which they capture the residual nonlinear trend.   In our first simulation, there is no nonlinear trend in the fitted values, and the first derivative is flat, so any curvature found in the fitted values will show as error in the second component of the RMSE.

Results from this decomposition can be found in Figure~\ref{fig:fd_linear} and Figure~\ref{fig:fd_nonlin} . Figure~\ref{fig:fd_linear} presents results for the low-dimensional linear trend while Figure~\ref{fig:fd_nonlin} presents results for residual non-linear trend terms.  We find that the proposed method performs well for both sets of effects. Interestingly in Figure~\ref{fig:fd_nonlin} we see that POLYMARS does more poorly in the non-linear setting but better in the discontinuous setting.\footnote{Beyond individual effects, sample average effects though, may be of interest to the researcher or policy maker.  For that reason, we also considered the level of the bias and generally find similar results though the Bayesian Bridge performed better than in the previous results.}

\begin{figure}[t!]
%\vspace{-.9in}
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{RMSE_FirstDerivative_linear.pdf}
    \vspace{-.25in}
  \end{center}\spacingset{1}
  \caption{{\bf A Comparison of RMSE Across Methods. Linear terms.}      } \label{fig:fd_linear}
\end{figure}

\begin{figure}[t!]
%\vspace{-.9in}
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{RMSE_FirstDerivative_nonlinear.pdf}
    \vspace{-.25in}
  \end{center}\spacingset{1}
  \caption{{\bf A Comparison of RMSE Across Methods. Residual non-linear terms}      } \label{fig:fd_nonlin}
\end{figure}



\subsection{Extrapolation Ranges and Expected Effects}
\label{sec:extrap}

We have so far focused on the partial derivative of an outcome with respect to a treatment as our primary estimand.  The researcher may want to know the causal effect of a larger shift in the treatment, or perhaps expected outcomes over a range of treatment values.  As a model of the conditional mean, MDE certainly {\it can}  estimate the outcome under any treatment level.  Estimates are only reliable, though, over the range of treatment values that have non-zero probability, $\supp(\widetilde T_i|X_i)$.  Doing so requires an estimate of the density of $\widetilde T_i|X_i$

We suggest using a bootstrap estimate.  Denote as $\widehat T_i$ the estimate of the treatment from our regression model;  $\widehat \epsilon_i$ the estimated residual; $\E(T_i|X_i)$ the population value;  $\widehat T_{ib}, \widehat \epsilon_{ib}$  estimated fitted values and residuals from a bootstrap replicate;  and $\widehat T_B=\{\widehat T_{ib}\}_{b=1}^B, \widehat \epsilon_B=\{\widehat \epsilon_{ib}\}_{b=1}^B$, and $q(A, \alpha)$ the $\alpha^{th}$ quantile of set $A$.

We cannot construct a bootstrapped estimate from $\widehat T_{ib} + \widehat \epsilon_{ib}$ since this value is $T_i$ for each replicate.  Instead, we construct our estimates from $T_{ib}=\widehat T_{ib} + \widehat \epsilon_{ib^\prime}$, for some $b \ne b^\prime$, breaking the correlation between fitted value and residuals for a given observation over bootstrap replicates.  Under the assumption that the error is independent of the systematic component, the set $T_B=\{T_{ib}\}_{b=1}^B$ can be used to construct a confidence interval.  For example, we can take use the percentile bootstrap, $[q(T_b, \alpha/2),q(T_B, 1-\alpha/2)]$ for false positive rate $\alpha$.  For proof, see Appendix \ref{app:bootstrap}.



\begin{comment}
A valid confidence interval can be constructed from $T_{ib}$.  To see this, assume $1-\alpha_T\%$ and $1-\alpha_\epsilon\%$ bootstrapped confidence intervals on $\widehat T_i$, $\widehat \epsilon_i$:
\begin{align}
\Pr(\E(\widetilde T_i|X_i) \in [T_l, T_u]) \ge 1-\alpha_T\\
\Pr(\widetilde T_i-\E(\widetilde T_i|X_i)\in [\epsilon_l, \epsilon_u]) \ge 1-\alpha_\epsilon.
\end{align}
Then,
\begin{align}
&\Pr(\widetilde T_i \in [T_l+\epsilon_l, T_u + \epsilon_u]) =\\
&\Pr(\widetilde T_i \in [T_l+\epsilon_l, T_u + \epsilon_u]|\widetilde T_i-\E(\widetilde T_i|X_i)\in [\epsilon_l, \epsilon_u] ) \Pr(\widetilde T_i-\E(\widetilde T_i|X_i)\in [\epsilon_l, \epsilon_u] )\ge\\
&\Pr(\widetilde T_i \in [T_l, T_u ] ) \Pr(\widetilde T_i-\E(\widetilde T_i|X_i)\in [\epsilon_l, \epsilon_u] )\ge
&(1-\alpha_T)\times (1-\alpha_\epsilon).
\end{align}
This shows that a $1-\alpha\%$ bootstrapped confidence interval

of $\E(T_i|X_i)$, the set of $B$ bootstrap estimates $\widehat T_B=\{\widehat T_{ib}\}_{b=1}^B$, and $q(\widehat T_B, \alpha)$ the $\alpha^{th}$ quantile of the bootstrap estimates.


Denote as $\widehat T_i$ the estimate of the treatment from our regression model,  $\E(T_i|X_i)$ the population value,  $\widehat T_{ib}$ a bootstrap estimate of $\E(T_i|X_i)$, the set of $B$ bootstrap estimates $\widehat T_B=\{\widehat T_{ib}\}_{b=1}^B$, and $q(\widehat T_B, \alpha)$ the $\alpha^{th}$ quantile of the bootstrap estimates.

We  use $\widehat T_B$ to construct an approximation to the distribution of $\widetilde T_i$, as $\widehat T_{ib}\ -\widehat T_{i} \convd \widehat T_{i}\ -\E(T_i|X_i)$.

This suggests a $1-\alpha\%$ symmetric extrapolation range as
\begin{align*}
\l[2\widehat T_i -q\l(\widehat T_B, 1-\frac{\alpha}{2}\r), 2\widehat T_i -q\l(\widehat T_B, \frac{\alpha}{2}\r) \r]
\end{align*}
with $\alpha$ some pre-specified level, say $0.05$ or $0.1$.
\end{comment}


This bootstrap estimate of the sampling distribution offers three separate advances.  First, we can use it to estimate a reasonable range of extrapolation for each covariate profile.  Doing so helps give a sense of what range of counterfactual predictions are supported by the data.

Second, the bootstrap estimated support can be used {\it ex ante}, say in an experiment, or {\it ex post}, in an observational study, to find observations with common support.  In an experimental study, a randomization with poor in-sample balance on  prognostic variable, can be rerun \citep[e.g.][]{Morgan12}.  In  a field study, where experimental units may be quite heterogeneous, the researcher may have concerns over the ability to compare the treated and control group in a given randomization.  The support for the treated and control can be estimated, and consequently compared, as a means of assessing whether the treatment and control groups are indeed comparable.  In an observational study, the estimated support can be used for trimming.  Rather than trim off the marginal covariates \citep[e.g.][]{king2006dangers}, the bootstrap estimates can be used for trimming off the actual estimated support of the treatment variable.

Lastly, the bootstrap estimate can be used to extend from a local effect at the observed value, $\lim_{\delta \rightarrow 0} \nabla_{T_i}(\delta)$, to an expected effect.  The local effect conditions on a particular value of the treatment, $\widetilde T_i=T_i$.  It could very well be, though, that the effect at the observed $T_i$ may be different than an expected value over all possible $T_i$.  The researcher interested in the expected effect over the treatment, $\E_{\widetilde T_i}(\lim_{\delta \rightarrow 0} \nabla_{T_i}(\delta) )$, can use $T_B$ to approximate this integral,
\begin{align}
\widehat \E_{\widetilde T_i}\l\{\lim_{\delta \rightarrow 0} \nabla_{T_i}(\delta)  \r\}& = \frac{1}{B} \sum_{b=1}^B   \lim_{\delta \rightarrow 0} \nabla_{T_{ib}}(\delta),
\end{align}
the expected effect of the treatment.


\section{Extensions}
\label{sec:extensions}

We next consider several extensions of the proposed method.  We first extend our local effect estimand to two common structural models: instrumental variable (IV) analyses and mediation analysis.  We show that MDE provides new insights into the problems.  In the IV setting, MDE lets us relax the monotonicity assumption, that the instrument uniformly encourage or discourage the treatment across the sample.  In the mediation setting, we show that MDE allows us to maintain a useful identity, that the total effect is additive in the direct and mediation effects, even with nonlinear models.  Next, we discuss a means for estimating the range of values for the treatment over which the model could predict.  We then use this approach to change our estimand to the expected value of the local effect over the range of treatment values, $\E_{\widetilde T_i}(\lim_{\delta \rightarrow 0} \nabla_{T_i}(\delta) )$.   Last, we show that MDE easily accommodates binary treatments, the most studied case in the literature.

\subsection{Instrumental Variables}

A treatment and outcome may be mutually determined, and thereby a correlative or regression analysis may not recover the causal effect of the treatment.  In these cases, an instrument, or ``encouragement,'' can  recover a causal effect.  The instrument is assumed to have a direct effect on the treatment but no direct effect on the treatment, thereby providing an experimental handle in the study.

Two  problems emerge:  the instrument may only affect some observations, the compliers, and not others; and the instrument may have a positive impact on some observations and negative on others.  These issues pose problems to both internal and external validity.  It is unclear which observations are actually impacted by the instrument and hence driving the causal effect estimate.  The second concern, that there are ``no-defiers,'' is assumed away in the binary treatment/instrument case, while this assumption is embedded in a linear first-stage specification.

To formalize, we assume  the treatment is not conditionally independent of the outcome given the covariates, $Y_i(\widetilde T_i) \notindep \widetilde T_i|X_i$, thereby biasing the  estimate of the causal effect.  We assume the existence of a pre-treatment instrument $\widetilde Z_i$ that helps resolve the issue.  The instrument, which follows law $F^Z_{X_i}$ and has observed value $Z_i$, enters  the potential outcome function for the treatment as $\widetilde T_i=T_i(\widetilde Z_i)$. For identification, we make the exclusion restriction that $Z_i$ has no direct effect on the outcome, $Y_i(T_i(\widetilde Z_i), \widetilde Z_i) = Y_i(T_i(\widetilde Z_i), Z_i^\prime)$ for all $Z_i^\prime \in \supp(\widetilde Z_i)$, or, equivalently, $Y_i(T_i(\widetilde Z_i), \widetilde Z_i)  \indep \widetilde Z_i|\widetilde T_i, X_i$.  The observed data is $[Y_i, T_i, Z_i, X_i^\top]^\top$.

In order to recover a consistent estimate of the causal effect of the treatment on the outcome, we  trace the exogenous variation of the instrument through the treatment and onto the outcome.  The {\it encouragement} from perturbing the instrument on observation $i$ is
\begin{align}
\nabla_{Z_i}^{IV}(\delta) = \E(T_i(Z_i+\delta)- T_i(Z_i) |X_i)\label{eq:encouragement},
\end{align}
where as above we have to condition on $X_i$ to recover an estimate.
The {\it intent to treat on the treated}, or {\it ITT}, is the impact of this encouragement on the outcome,
\begin{align}
\nabla_{T_i}^{IV}(\delta)=  \E\l\{Y_i\l(T_i(Z_i)+\nabla_{Z_i}^{IV}(\delta)\r) -Y_i\l(T_i(Z_i)\r)|X_i \r\}\label{eq:itt}.
\end{align}

The effect $\lim_{\delta \rightarrow 0} \nabla_{T_i}(\delta)=\E(d Y_i(T_i)/d T_i|X_i)$ is not identified because of the failure of ignorability.   Instead target the causal quantity
\begin{align}
&\E\l(\frac{\partial Y_i(T_i(Z),Z)}{\partial Z} \bigg|X_i,Z= Z_i \r)=\\
&\E\l(\frac{\partial Y_i(T_i(Z))}{\partial Z} \bigg|X_i,Z= Z_i \r)=\\
&\E\l(\frac{\partial Y_i(T_i)}{\partial T_i} \frac{\partial T_i(Z)}{\partial Z} \bigg|X_i,Z= Z_i \r)
\end{align}
where the second line comes from the exclusion restriction and the third by the chain rule.  We then use the continuous mapping theorem to give
\begin{align}
\E\l(\frac{\partial Y_i(T_i)}{\partial T_i}\bigg|X_i, Z=Z_i\r)&=\plim \frac{\E(\partial Y_i(T_i(Z))/\partial Z|X_i,Z= Z_i )}{\E(\partial T_i(Z)/\partial Z\big|X_i,Z= Z_i )}\\
&=\lim_{\delta \rightarrow 0 }  \frac{ \nabla_{T_i}^{IV}(\delta)}{ \nabla_{Z_i}^{IV}(\delta)}\label{eq:lice},
\end{align}
which we refer to as the Local Instrumental Causal Effect (LICE) and denote $\lim_{\delta \rightarrow 0} \nabla_{T_i,Z_i}^{IV}(\delta)$. Like with TSLS, the plug-in estimate is a biased but consistent estimate of the desired causal effect.



%\paragraph{Comparison with Two-Stage Least Squares Estimate.}

The LICE differs in an important way from the standard Wald or two-stage least squares estimates (TSLS).  The TSLS estimate under a linear structural equation model, which reduces to the Wald estimate in the binary treatment/encouragement seeting,  is calculated from sample covariances as
\begin{align}
\widehat \theta^{TSLS} = \frac{\widehat \Cov_{\mathcal S}(Y_i, Z_i|X_i)}{\widehat \Cov_{\mathcal S}(T_i, Z_i|X_i)},
\end{align}
a ratio of sample covariances, after $X_i$ has been partialed out. The TSLS estimator equals the OLS estimate of the effect of $T_i$ on $Y_i$ when $Z_i=T_i \;\forall\; i$, i.e. if encouragement is perfect.  Compared to our estimate in Equation \ref{eq:lice}, we see that the TSLS estimate averages the numerator and denominator over compliers, defiers, and observations for which no causal effect is identified.

The estimated LICE is a TSLS estimate, but one that varies with $X_i$, thereby  bringing the underlying heterogeneity to the fore.  To see this, consider
\begin{align}
\lim_{\delta \rightarrow 0} \nabla_{T_i,Z_i}^{IV}(\delta) &= \lim_{\delta \rightarrow 0}  \frac{\nabla_{T_i}^{IV}(\delta)}{\nabla_{Z_i}^{IV}(\delta)}\\
&=\frac{\E(\partial Y_i(T_i(Z))/\partial Z|X_i,Z= Z_i )}{\E(\partial T_i(Z)/\partial Z\big|X_i,Z= Z_i )}\\
&=\frac{\Cov(Y_i, Z_i|X_i)/\Var(Z_i|X_i)}{\Cov(T_i,Z_i|X_i)/\Var(Z_i|X_i)}\\
&=\frac{\Cov(Y_i, Z_i|X_i)}{\Cov(T_i,Z_i|X_i)},
\end{align}
where going from the second to the third line follows from the exclusion restriction and we assume that the outcome and treatment functions are differentiable in the instrument.

%\footnote{This is closely related to, yet differs from, \citet{Heckman07a, Heckman07b, Heckman05, Heckman99}. We estimate a { \it pointwise} effect rather than integrating over a first-stage choice set.  Also, we are not limited to the binary instrument/treatment.  See \citep{kennedyiv} for concurrent, and fascinating, work. }





\begin{comment}
Our desired causal estimate is
\begin{align}
which we estimate using
\begin{align}
\frac{\widehat \E\l(\lim_{\delta \rightarrow 0 } \nabla_{T_i}^{IV}(\delta)|X_i\r)}{\widehat\E\l( \lim_{\delta \rightarrow 0 } \nabla_{Z_i}^{IV}(\delta)|X_i\r)}\label{eq:lice2},
\end{align}
which, like a TSLS estimate, is a biased but consistent estimate of the desired causal effect.

\E\l(\frac{d}{dT}Y_i(T_i)|X_i\r)=\\
\E\l(\frac{\partial Y_i(T_i(Z_i),Z_i)/\partial Z_i}{\partial T_i(Z_i)/\partial Z_i}|X_i \r)=\\
\E\l(\frac{\partial Y_i(T_i(Z_i))/\partial Z_i}{\partial T_i(Z_i)/\partial Z_i}|X_i \r)
\end{align}
\end{comment}


   Since our estimand is conditioned by $X_i$, we do not need to assume homogeneity in the denominator (monotonicity) across values of $X_i$ in order to estimate into which principal stratum each observation falls  \citep[as required by][]{Abadie2003, aronow2013beyond}.  We will refer to observations with a positively, negatively, and zero impact of the instrument on the treatment as compliers, defiers, and non-compliers, respectively \citep{angrist1996identification}.   Were we to assume latent index models underlying each  a binary treatment and binary instrument, we would recover the same model used in the Local IV (LIV) estimate in \citet{Heckman07a, Heckman07b, Heckman05, Heckman99}.  MDE is  means of estimation, and as such, either the LIV or LICE can be estimated; we defer a discussion comparing the two estimands to below.

For estimation, we  fit two models: one for the treatment as a function of the instrument and a second for the outcome as a function of the treatment.  We do so using the four-step procedure outlined above, and then use the plug-in estimates for equations \ref{eq:encouragement} and \ref{eq:itt} to estimate the LICE in equation \ref{eq:lice}.

\paragraph{Threshold for Estimating Compliers}  The LICE only exists for observations for which the encouragement is non-zero.  To estimate these observations, we implement a feasible estimate of the  Oracle bound, as given in Inequality \ref{eq:oracle2}, to allow us to differentiate observations with some systematic perturbation from encouragement estimates that are likely noise.
Since  the Oracle bound does not vary across the sample, we implement an adaptive bound that narrows where we estimate signal and expands in a noisy region, constructed from a pointwise estimated degree of freedom.

\citet{Stein81} showed in the normal regression model $Y_i \sim \mathcal N(\widehat \mu_i(Y_i),  \sigma^2)$ for generic conditional mean function $\widehat \mu_i(Y): Y \rightarrow \widehat Y_i$, a degree of freedom estimate can be recovered as
\begin{align}
\widehat {df}_i=\frac{\partial \widehat \mu_i(Y_i)}{\partial Y_i}=\frac{\Cov(Y_i, \widehat \mu_i(Y_i))}{\sigma^2}
\end{align}
and the model degree of freedom can be estimated as $\widehat{df}=\sum_{i=1}^N \widehat{df}_i$.  This definition coincides with the trace of the projection matrix when $\widehat \mu_i(Y_i)$ is linear in $Y_i$.



Decomposing into parts of the design that covary with $T_i$ and those that do not, we can decompose the degrees of freedom into
\begin{align}
\widehat {df}_i&=\frac{\partial [R^X(X_i): R^{TX}(T_i, X_i)]^\top \widehat c}{\partial Y_i}\\
&=\frac{\partial [R^X(X_i)]^\top \widehat c^X}{\partial Y_i}+\frac{\partial [R^{TX}(T_i, X_i)]^\top \widehat c^{TX}}{\partial Y_i}.
\end{align}
We then take the degree of freedom estimate associated with $\widehat \nabla_{T_i} (\delta)$ as
\begin{align}
\widehat {df}^\nabla_i = \frac{\partial [R^{TX}(T_i, X_i)]^\top \widehat c^{TX}}{\partial Y_i}.
\end{align}
This estimate will be larger the more signal there is at each observation.

The adaptive component of our threshold is of the form
\begin{align}
\widehat df^{adapt}_i=\frac{1/n}{1/n+\widehat {df}^\nabla_i }
\end{align}
which takes on a value of $1$ when there is no signal for observation $i, (\widehat {df}^\nabla_i =0)$.  If the true model is additive and linear in the treatment, we would  see $\widehat {df}^\nabla_i =1/n$, which gives $\widehat df_i^{adapt}=1/2$.  As the model grows more complex at a given point, the threshold will shrink.

We combine the degree of freedom bound and oracle estimate in our threshold as
\begin{align}
\widehat{ \mathbf 1}\l(\nabla_{Z_i}^{IV}(\delta_i^{IV}) \ne 0\r)=\mathbf 1\l(\big|\widehat \nabla_{Z_i}^{IV}(\delta_i^{IV}) \big|>
C \frac{\widehat \lambda ||\widehat w_k||_\infty \widehat \sigma}{n \widehat \phi_0}  \widehat df_i^{adapt} \r)
\end{align}
with $||\widehat w_k||_\infty$ the largest weight.  $C$ is a user-selected constant, but we show in our simulation that taking $C=2$ helps select observations impacted by the instrument. Lastly, the compatibility constant is the smallest eigenvalue of the covariance of the true model predicting the gradient.  We estimate it using columns of the submatrix of the design that contains interactions with the treatment.  We find this matrix may be ill-conditioned or even rank-deficient, so we estimate the smallest eigenvalue such that the eigenvalues above it explain $90\%$ of the variance in the selected model.

The geometry of this threshold incorporates a ``complexity penalty'' in the selection process.  If the estimated model for the treatment effect is simple, the design will be low-dimension with a flat spectrum, say a linear model.  As the model grows more complex, the design will incorporate more nonparametric bases and its smallest eigenvalues will be closer to zero, inflating the threshold.  The adaptive component, $\widehat df_i^{adapt}$, serves to adjust for local effects.

We include simulation evidence for the proposed threshold's efficacy and reliability in the supplemental materials (Appendix~\ref{app:ivsims}).


\subsection{Causal Mediation}

Mediation analysis examines the pathway through which a treatment variable impacts some outcome through an intermediate variable, a mediator. Mediation analysis hence examines mechanisms rather than treatment effects; for a review, see  \citet[][]{vanderweele2015explanation}.

A key insight of mediation analysis is that the total treatment effect decomposes into the direct effect plus the mediated effect.  In a linear SEM framework,  the mediated effect is the product of the mediator direct effect and intermediate direct effect \citep{mackinnon2002comparison}, and it can be estimated as a product of coefficients. \citet{imai2010general} show that the product rule  does not apply to estimated coefficients unless the underlying model is linear.  We show that MDE allows us to maintain this key identity, as we are fitting a model that is a local linearization around the treatment and mediator at each point.

Consider first a post-treatment mediator equipped with its own potential outcome function, $\widetilde M_i=M_i(\widetilde T_i)$ with law $F^M_{T_i, X_i}$ and observed value $M_i = M_i(T_i)$.  The observed outcome is now $Y_i = Y_i\l(M_i, T_i\r)$ and the observed data is $\l[Y_i, M_i, T_i, X_i^\top\r]^\top$.  Our goal is to differentiate three effects: the total effect of the treatment on the outcome, the direct effect of the treatment on the outcome, and the mediated effect of the treatment through the mediator on the outcome.  For identification, we will assume that the covariates are sufficiently rich to render the outcome, mediator, and treatment sequentially ignorable \citep{imai2010identification}:  $Y_i(\widetilde M_i, \widetilde T_i) \indep M_i(\widetilde T_i)  |\widetilde T_i, X_i$ and $\l\{Y_i\l(\widetilde M_i, \widetilde T_i\r), M_i(\widetilde T_i) \r\}\indep \widetilde T_i|X_i$.  Second, we assume that all considered manipulations of mediator and treatment are on the support of $\widetilde M_i|X_i$ and $\widetilde T_i|X_i$.


We then estimate the observation-level total effect, direct effect, and mediated effect as
\begin{align}
\textrm{Total effect:}\;\; \nabla_{T_i}^{TE}(\delta_i^{TE}) &=\frac{1}{\delta_i^{TE}}\E \l\{Y_i\l( M_i(T_i+\delta_i^{TE}), T_i+\delta_i^{TE} \r)-Y_i\l(M_i, T_i \r)|X_i\r\}\label{eq:eff1}\\
\textrm{Direct effect:}\;\;\nabla_{T_i}^{DE}(\delta_i^{DE}) &=\frac{1}{\delta_i^{DE}} \E\l\{Y_i\l( M_i, T_i+\delta_i^{DE} \r)-Y_i\l(M_i, T_i \r)|X_i\r\}\label{eq:eff2}\\
\textrm{Mediated effect:}\;\;\nabla_{T_i}^{ME}(\delta_i^{ME}) &=\frac{1}{\delta_i^{ME}} \E\l\{Y_i\l( M_i(T_i+\delta_i^{ME}), T_i\r)-Y_i\l(M_i, T_i \r)|X_i \r\}.\label{eq:eff3}
\end{align}
  We will also make use of
\begin{align}
\textrm{Mediator direct effect:}\;\;\nabla_{M_i}(\delta_i^{M}) &=\frac{1}{\delta_i^{M}}\E \l\{Y_i\l(M_i+\delta_i^M, T_i\r)-Y_i\l( M_i,T_i \r)|X_i\r\}\label{eq:eff4}\\
\textrm{First stage mediation effect:}\;\;\nabla_{T_i}(\delta_i^{T}) &=\frac{1}{\delta_i^{T}} \E\l\{M_i\l( T_i+\delta_i^{T}\r)-M_i\l( T_i \r)|X_i\r\}\label{eq:eff5}
\end{align}

If we assume that the potential outcome functions $Y_i(\cdot)$ and $M_i(\cdot)$ are differentiable in their manipulable arguments, we see that the total effect decomposes into a sum of the direct and mediated effects. The law of the total derivative gives
\begin{align}
\lim_{\delta_i^{TE}\rightarrow 0} \nabla_{T_i}^{TE}(\delta_i^{TE})&= \lim_{\delta_i^{DE}\rightarrow 0}\nabla_{T_i}^{DE}(\delta_i^{DE}) + \lim_{\delta_i^{M}\rightarrow 0}\nabla_{M_i}(\delta_i^{M})\times \lim_{\delta_i^{T}\rightarrow 0}\ \nabla_{T_i}(\delta_i^{T}).
\end{align}
The second summand on the righthand side can be simplified by the chain rule as
\begin{align}
\lim_{\delta_i^{TE}\rightarrow 0} \nabla_{T_i}^{TE}(\delta_i^{TE}) &= \lim_{\delta_i^{DE}\rightarrow 0}\nabla_{T_i}^{DE}(\delta_i^{DE}) + \lim_{\delta_i^{ME}\rightarrow 0}\ \nabla_{T_i}^{ME}(\delta_i^{ME}).
\end{align}
This gives us two well-known results: the mediated effect is the product of the mediator direct effect and intermediate direct effect \citep{mackinnon2002comparison} and that the total effect is additive in the direct and mediated effects.

In our framework we thus estimate two models, one for the outcome and one for the mediator.  The models that we fit are
\begin{align}
Y_i &= \mu_Y+R_Y(M_i, T_i, X_i)^\top c_Y + \epsilon_i^Y\\
M_i &= \mu_M+R_M( T_i, X_i)^\top c_M + \epsilon_i^M
\end{align}
where, again, the $R_\cdot(\cdot)$ vectors are a post-screened set of bases.  We can now recover estimates of the effects in  Equations \ref{eq:eff1}--\ref{eq:eff5}.

Interestingly, while \citet{imai2010general} show that the product rule above does not apply to estimated coefficients unless the underlying model is linear, we obtain a different result. Since we are using a local linearization around each point as a function of the treatment and mediator, the product rule can be extended.


\subsection{Controlled Direct Effects and Sequential g-estimation}



A natural extension of our method is to the case where controlled direct effects are of interest and there are intermediate confounders \citep[as in][]{acharya2015explaining,vansteelandt2009estimating}. We focus on estimating a controlled direct effect, but show how MDE can then generalize to a far wider class of structural mean models.

The controlled direct effect is the effect of a shift in treatment holding the mediating variable fixed at some level $m$:
\begin{align}
\textrm{Controlled direct effect:}\;\; \nabla_{T_i}^{CDE}(\delta_i^{CDE}) &=\frac{1}{\delta_i^{CDE}}\E \l\{Y_i\l( M_i=m, T_i+\delta_i^{CDE} \r)-Y_i\l(M_i=m, T_i \r)|X_i\r\}
\end{align}
Identification assumptions are nearly identical to the mediation case above, except the CDE allows for post-treatment confounders in the condition set.  Decompose $X_i=[X_i^{pre}:X_{i}^{post}]$ for whether the covariates are observed pre- or post-treatment. The assumptions are, first, $Y_i(\widetilde M_i, \widetilde T_i) \indep M_i(\widetilde T_i)  |\widetilde T_i, X_i^{pre}, X_i^{post}$ and $\l\{Y_i\l(\widetilde M_i, \widetilde T_i\r), M_i(\widetilde T_i) \r\}\indep \widetilde T_i|X_i^{pre}$, and, second, all manipulations are on the support of $\widetilde M_i|X_i^{pre},X_i^{post}$ and $\widetilde T_i|X_i^{pre}$.

The estimate of the controlled direct effect is then simply
\begin{align}
\widehat \nabla_{T_i}^{CDE}(\delta_i^{CDE}) &=\frac{1}{\delta_i^{CDE}}(R_i(M_i=m,T_i+\delta_i^{CDE},X_i)-R_i(M_i=m,T_i+\delta_i^{CDE},X_i))^\top \widehat c
\end{align}

We note two advances offered by our estimation strategy.  First, estimation can be done in one step since we need not model the mediator as a function of treatment, as we did in recovering the mediation effect above.  The method of sequential  $g$-estimation, most commonly used in these models, is a multi-stage method, where summaries from earlier models are used as a bias-correction in later models.  For MDE, the bias correction is incorporated through the bases: estimation utilizes bases $R_i(M_i, T_i,X_i)$, but prediction at $R_i(m, T_i, X_i)$.  Second, we can relax the assumption that post-treatment confounders are independent of the mediator, which is required for the simplest form of $g$-estimation.

More generally, structural mean models can be estimated using the basic framework of MDE.  The basic element of these models, the ``blip function,'' is the conditional counterfactual difference in outcome between two observations.  Under an identity link, which we consider in this paper, the blip function for a point-treatment is
\begin{align}
\psi^\ast\l((T_i^{\prime\prime},T_i^{\prime}),X_i;c^o\r)=\E(Y_i(T_i^{\prime\prime})-Y_i(T_i^{\prime})|X_i)
\end{align}
where $c^o$ parameterizes the true mean function.  This blip function is constructed to satisfy
\begin{align}
\E\l\{Y_i(T_i^{\prime\prime}) -\psi^\ast\l((T_i^{\prime\prime},T_i^{\prime}),X_i;c^o\r)|X_i\r\}=\E\l\{Y_i(T_i^{\prime})|X_i\r\}
\end{align}
allowing recovery of an estimate of the outcome at some desired value, $T_i^\prime$.  This technology can be extended to recover consistent treatment effect estimates across a variety of structural settings \citep[see][ for an excellent review]{vansteelandt2009estimating}.

In our notation, this blip function $\psi^\ast\l((T_i^{\prime\prime},T_i^{\prime}),X_i;c^o\r)=\E\l\{\Delta_i\l(T_i^{\prime\prime},T_i^{\prime}\r))|X_i\r\}$, a generalization of our primary parameter of interest, $\nabla_{T_i}(\delta)$.  The methods described above show how MDE can recover a blip function but also assess whether the desired conditional in the blip down has positive support, a particular concern with a continuous treatment where structural mean models and their extensions are most useful \citep[][p. 23]{vansteelandt2009estimating}.

\begin{comment}
Estimation in this model is actually simpler than the estimating the natural direct effect above because we need not consider model the mediator as a function of treatment.  Modeling the

MDE extends to additional

In these contexts a two step sequential g-estimation procedure is used. Central to estimating controlled-direct effects, and working with structural mean models in general, is the ``blip function,''

The first step regresses the outcome on the treatment, mediator, and all covariates including those that are pre-treatment and intermediate. The model could be extended to included interactions between the treatment and any pre-treatment covariates. Next, a de-mediated function (or ``blip down'') function is calculated from this equation. This is just the coefficient on the mediator, multiplied by each observation's value of the mediator, and subtracted from the outcome. The second step then regresses the demediated outcome on the treatment and pre-treatment confounders, and the coefficient on the treatment is the estimate of the controlled direct effect.

\citet{acharya2015explaining} note the strong modelling assumptions required to use this method when dealing with continuous variables (non-parametric identification is shown for binary or limited value treatment and mediators variables in \citet{robins2000marginal}). We can relax these modeling assumptions, such as linearity and additivity, by using our nonparametric regression to ``de-mediate'' the outcome.  With the de-mediated outcome, we could then use the method of direct estimation to for estimating the impact of a controlled direct fluctuation of the treatment on the outcome. With this connection to controlled direct effect we open up a broader connection to the literature on sequential g-estimation.

\end{comment}


\subsection{Bagging}
We have several reasons to be concerned that our estimates may have a large sampling variance.  The LASSO problem minimized in Equation \ref{eq:lassoderiv} conditions on the estimated residuals.  This suggests that results may be sensitive to sample-specific outlying residuals.  Second, we are fitting a high-dimensional nonparametric regression, which is known to generate high-variance estimates near the boundaries of the covariate space.   Beyond the generic issues of nonparametric regression, several of our examples and simulations involve cases where we confront highly-skewed or fat-tailed distributions.  For example, in an instrumental variable setting, the effect estimate is a ratio estimator that may have such fat tails that the sampling distribution may have no finite moments.  In one of our empirical examples, the outcome is earnings, which is highly right-skewed.


To reduce the variance of our estimates, we turn to bootstrap aggregation \citep[bagged,][]{Buhlmann02}.  We implement bagging use a Rademacher-wild bootstrap, to account for possible heteroskedasticity \citep[e.g.,][]{Davidson08}.  Bagging normally involves taking the mean across bootstrapped samples, but as we are worried about skew or whether the mean is even finite, we take as our bagged estimate the median across bootstrap samples for each observation.  We show below that the method can lead to a decrease root-mean-squared error with only little increase in bias.

We confront a similar problem with our IV estimator.  The LICE is a ratio estimator, and we worry that the estimator's sampling distribution may be Cauchy or approximately so.  In order to stabilize the estimation, we utilize median bagged estimates,
\begin{align}
\begin{split}
\widehat \nabla_{Z_i}^{IV;boot}(\delta_i^{IV}) =\med_B\l( \frac{1}{\l\{R_T(Z_i+\delta_i^{IV}, X_i) - R_T(Z_i,X_i)\r\}^\top \widehat c_{T,b}}\r)\times \\\med_B\l(\l\{R\l(R_T(Z_i+\delta_i^{IV}, X_i) ^\top \widehat c_{T,b}, X_i\r) - R\l(R_T(Z_i, X_i) ^\top \widehat c_{T,b}, X_i\r) \r\}^\top\widehat c_b\r)
\end{split}
\end{align}
where $\med_B(a)$ refers to the median over bootstrap samples $b \in \{1, 2, \ldots, B\}$.  This median-bagging is effective at reducing the impact of wild or erratic estimates.





\subsection{Binary Treatments}

We have so far focused on and emphasized the method's utility with a continuous treatment variable, but there is nothing in our framework that prevents us from predicting counterfactuals with a binary or categorical treatment variable.  We next extend MDE to the case of a binary treatment, where most work in the causal literature has taken place, and show that it provides a reasonable means for effect estimation even relative to methods designed explicitly for the binary treatment case.

In our notation, the individual causal effect under a binary treatment can be written as
\begin{align}
\E(Y_i(1)-Y_i(0)|X_i)=\nabla_{T_i}(1-2\times T_i)=\E(\Delta_i(1,0)|X_i)
\end{align}
The most common estimands with a binary treatment are the sampe average treatment effect (ATE) and sample average treatment effect on the treated (ATT),
\begin{align}
ATE=\frac{1}{n} \sum_{i=1}^n \nabla_{T_i}(1-2\times T_i) = \\
ATT=\frac{1}{n_T} \sum_{i=1}^n \l(\nabla_{T_i}(1-2\times T_i)\r)\times \mathbf 1(T_i=1)
\end{align}
where $n_T= \sum_{i=1}^n  \mathbf 1(T_i=1)$.  Estimation can now proceed as described above, except we now consider differences instead of partial derivatives.



\section{Relationship to Earlier Work}
\label{sec:earlier}

Our approach to modeling counterfactuals incorporates insights from several different fields, as we discuss next.

\paragraph{Relationship to causal inference} Causal estimation generally involves a two-step  procedure.  In the first, a model of the treatment is used to characterize any confounding.  In the second, some summary statistic from the first stage, such as a the density estimate or conditional mean of the treatment, is incorporated into the outcome model through matching, subclassification, or inverse weighting.  See \citet{rubin1974estimating, rose:rubi:83, rose:rubi:84, robi:rotn:zhao:94, Pearl2000, tmlebook} for seminal work.

Misspecification of this treatment model is both inevitable and impactful \citep{kang2007demystifying}.   We sidestep the entire endeavor and directly target a fully saturated, nonparametric conditional mean function. In doing so we face efficiency reductions compared to a model that incorporates a correctly specified propensity score model \citep[e.g.,][]{robins1997toward}, yet this efficiency loss must be balanced against the modeling uncertainty that comes from not knowing the true model.  We show that focusing attention on the outcome model can provide a feasible, robust, and powerful means for causal estimation.

The bulk of the literature on causal inference has focused on the case of the case of binary or categorical treatment regimes \citep[For exceptions see ][]{imai2012causal,Hirano2005}.
We share some of the motivation in \citet{Austin2012, Hill2011, Lam13}, though these works focused on the binary treatment and did not extend the structural models.  We find below that the methods advocated by several of these works do not perform as well as our tensor regression model.  We also share a motivation with the targeted maximum-likelihood approach of \citet{tmlebook}.  The authors use a cross-validation tuned ensemble of methods, or ``Superlearner,'' to predict the treatment density, then incorporate these estimates in the second stage model so as to achieve semiparametrically efficient estimates of a treatment parameter.  Unlike this method, we target observation-level estimates rather than an aggregate parameter.  We also find the Superlearner performs poorly in estimating the derivative of the loss function, as minimizing the predictive risk may result in poor estimates of a derivative or treatment effect; see \citet{Athey16} and \citet{Horowitz14}, esp. sec 3.1.


\paragraph{Relationship to High-Dimensional Regression Modeling} Our spline basis is closest to the nonparametric specification in the ``POLYMARS'' multivariate tensor-spline model of \citet{Stone1997}. We differ from POLYMARS in two regards. First, our screening step allows us to include an ``ultra-high'' number of candidate bases \citep{Fan:2008}.  We  include spline bases of multiple degrees and interactions, then reduce the millions of possible bases to hundreds or thousands. Second, rather than conduct Rao and Wald tests for basis inclusion/exclusion, we use a sparse model to fit all of the screened bases at once \citep{RatkovicTingley17}.  Third, unlike existing sparse Bayesian models \citep[e.g.][]{Polson12, Rockova17}, we use a frequentist minmax argument to motivate our hyperprior selection as a function of $n,p$. We also utilize bagging in order to guard against erratic results attributable to skewed or fat-tailed sampling distributions \citep{Buhlmann02}.

We were inspired by work on tensor-product and smoothing spline ANOVA models \citet{Gu2002, Wood16, Wood06, Currie06, Eilers96}.  We differ from these methods primarily in scale.  We are considering tens or hundreds of covariates, a combinatorially growing number of tensor products, while focusing on the impact of only a single treatment.  Due to the complexity of the problem, we do not specify a penalty function or reproducing kernel, but instead construct tensor-product bases and let the sparse LASSO prior manage the regularization.  While this leads to inefficiencies in estimation, current software cannot handle the number of variables and tensor-interactions we are considering here.


\paragraph{Estimation of first derivatives.} In the case of a continuous or binary treatment, our method reduces to estimating the partial derivative or subderivative, respectively, of the response function with respect to the treatment variable. The econometric literature has long recognized this manipulation-based, and hence causal, interpretation of a structural equation model \citep{Haavelmo43, hanseneconometrics, Pearl2014, angr:pisc:09}.  We differ by targeting observation-level counterfactuals rather than a parameter in a structural model.  Even absent a causal interpretation,  estimating derivatives has been motivated by problems in engineering as well as economics; see \citet{OSullivan86,Horowitz14} for overviews and \citet{Hardle89, Newey94} for seminal work.  We differ primarily in estimation, through considering a large-$p$ setting as well as models that are nonparametric in all covariates. These methods also stay silent on how to choose a basis, whether $B$-splines, Hermite polynomials, radial basis functions, and so on.  We leverage our ignorability assumption, which is central to causal inference, to show that the appropriate basis function for a causal effect is locally linear almost everywhere, i.e. the piecewise linear cubic-spline or B-spline basis.



\paragraph{Relation to nonparametric and high-dimensional structural models}  Recent work has extended nonparametric and high-dimensional estimation to the instrumental variables setting, through either series estimation or LASSO selection \citep[][]{newey2013nonparametric, belloni2012sparse, Newey03}.  Like these methods, we use a regularized series estimator to recover a conditional mean.  We differ from these methods in that we construct our estimator from observation-level predictions under counterfactual manipulations, rather than targeting structural parameters.  Our work is perhaps closest to \citet{Heckman07a, Heckman07b, Heckman05, Heckman99}\begin{footnote}{ See \citet{kennedyiv} for concurrent, and fascinating, work. }\end{footnote}, \citet{athey2016solving}, and \citet{taddydeep}.   Heckman and Vytlacil develop the Local IV (LIV) in the binary treatment/encouragement setting as  the impact of the treatment on those indifferent between control and treatment at a given level of encouragement.  The LICE considers the impact of a perturbation to the encouragement at the observed data.  In a policy setting, where the encouragement can be controlled, the LIV may be preferred, while in an observational setting, where interest is on the data as observed, the LICE is likely to be preferred.\footnote{See Section~\ref{sec:extrap} for one way to bridge the two estimands. One could for example calculate the set of observations for which a given level of encouragement could be supported.  We also note that, as with our approach the LIV can be aggregated to an average effect. For the Heckman and Vytlacil framework this follows precisely because it is a latent index model with an assumed error (a probit model).  In the continuous setting, like we consider here, that expectation requires either assuming the error density or estimating it, as we do.}  We emphasize, though, that MDE can be used to estimate either the LICE or LIV, using a rich mean model unavailable to these previous papers, and handles settings beyond a binary treatment/encouragement.




In additional related work,  \citet{athey2016solving} use a random forest to estimated a kernel density for each observation, then use the density weights in a two-stage least-squares calculation.  As with us,  \citet{athey2016solving} generate observation-level effect estimates and use a high-dimensional model to avoid the curse of dimensionality.  Like \citet{taddydeep}, we use a flexible regression model to estimate conditional mean functions.  While we generate pointwise counterfactuals at the two stages, \citet{taddydeep} uses the first stage for a second-stage residual correction, which is the ``control function'' approach described in \cite{Horowitz14}.   Our primary additions over  \citet{athey2016solving} and  \citet{taddydeep}  come from both returning a first-stage estimate, which offers important insight into internal validity, and using the oracle bound to identify non-compliers in the data.  The importance of separate estimates at each stage will allow plug-in estimates for other causal structural models \citep{vanderweele2015explanation}, as we illustrate with our discussion of mediation and sequential-g estimation above.   Lastly, we are not the first to note that the sampling distribution of instrumental variable estimates can be erratic and non-normal \citep{bound1995problems, imbens2005robust}.  Rather than turning to rank-based estimates, we instead implement a bagging procedure  to smooth over a possibly erratic sampling distribution.

%These works involve a behavioral foundation, and the authors advocate averaging over the entire choice set, whereas we focus on the local responsiveness at the observed treatment level.

\section{Analyses}
\label{sec:app}


\subsection{Application: Binary Treatment}

As a validation exercise, we test MDE against several existing methods using data  first put forward by \citet{lalonde1986evaluating}, but see also  \citet[][]{smith2005does, Diamond2012, CBPS14}.   The data consist of an experimental subset and and observational subset.  The experimental subset contains the results from a policy experiment, the National Supported Work Study, with participants randomly assigned to a treated ($n_T=297$) and untreated ($n_C=425$) group.  The treatment consisted of a job-training program and participants were hard-to-employ individuals across 15 sites in the United States in 1976.  The outcome is 1978 earnings and observed covariates include the participants' age, years of schooling, whether the individual received a high school degree, indicators for race (black, hispanic), an indicator for whether the participant is married, previous earnings from 1974 and 1975, and indicators for whether 1974 or 1975 earnings were zero.\begin{footnote}{The LaLonde data contain several subsets that have been used as a benchmark analysis.  One subset, analyzed by  \citet{dehejia1999causal} subsets on the outcome variable, returning a data on which most methods perform quite well.  We focus on the full original experimental data, which poses a greater challenge.  }\end{footnote}  The observational dataset comes from the Panel Study for Income Dynamics, a survey of low-income individuals, with data on the same covariates as the experimental data $(n_P=2915)$.

We conduct three separate tests. The first uses the experimental treated and observational untreated, with the goal of recovering the experimental result ($\$886.30$). This is perhaps the most policy-relevant comparison, as it may allow for methods that can estimate the causal impact of a policy intervention in the absence of a randomized control group.  The second test is a placebo test, considering the experimental and observational untreated groups.    As no one in the treated or untreated group received the treatment, the known true effect is zero, and any non-zero estimate has been termed {\it evaluation bias} \citep{smith2005does, CBPS14}.  In the third test, we compare the experimental treated to the experimental untreated, in order to assess the extent to which a causal estimation method can recognize experimental data and recover an effect close to the simple difference-in-means.


We compare the proposed method, in both its point estimate and bagged implementations, to several existing methods. We include the horseshoe model and cross-validated LASSO, both using our nonparametric bases.  We  also include BART, gradient boosting (GBM), and POLYMARS, an existing nonparametric spline model.    We include three methods designed for a binary treatment: logistic propensity score matching \citep[Propensity,][]{ho:imai:king:stua:07}, the covariate balancing propensity score \citep[CBPS,][]{CBPS14} and the targeted maximum likelihood estimate \citep{gruber2011tmle}.

% latex table generated in R 3.3.1 by xtable 1.8-2 package
% Wed Jan 11 20:17:13 2017
\begin{table}[tb!] \spacingset{1}
\centering
\begin{tabular}{r|rrr:r}
  \hline

Treatment Group &  Experimental & Experimental  & Experimental &Total \\
Control Grooup &Observational& Observational&Experimental &Absloute Bias\\
  \hline
  Truth & 886.30 & 0.00 & 886.30&0 \\
  \hdashline
  MDE & 850.82 & 190.55 & 572.01 & 540.31 \\
    MDE, Bagged & 960.16 & 96.66 & 416.90 & 639.92 \\
  CBPS & 453.19 & -300.28 & 872.27 & 747.42 \\
  TMLE& -128.11 & -520.48 & 739.26 & 1681.93 \\
  Horseshoe& 901.67 & 1600.79 & -93.01 & 2595.47 \\
 BART& -477.87 & -1265.07 & 803.77 & 2711.77 \\
  GBM & -864.43 & -1402.73 & 866.10 & 3173.66 \\
  POLYMARS & -228.30 & -1863.41 & 0.00 & 3864.32 \\
  LASSO & -1204.34 & -1726.13 & 944.94 & 3875.41 \\
  Propensity & -3728.93 & -5597.30 & 1067.96 & 10394.19 \\
   \hline
\end{tabular}\caption{{\bf Treatment Effect Estimates, LaLonde Data.}   The treated and control groups used in the comparison are given in the top two rows. Columns consist of the target, either the experimental benchmark $(\$886.30)$ or zero.  Results from each method are then listed below and the final column gives the absolute bias across comparisons.  Methods are listed in order of performance on this final measure.  We see that both MDE and the bagged variant perform well across the comparisons.  }\label{tab:lalonde1}
\end{table}


Results from this analysis can be found in Table \ref{tab:lalonde1}.  The treated and control groups used in the comparison are given in the top two rows. Columns consist of the target, either the experimental benchmark $(\$886.30)$ or zero.  Results from each method are then listed below and the final column gives the absolute bias across comparisons.  Methods are listed in order of performance on this final measure.  We see that MDE performs well relative to other methods.  CBPS achieves a bias across the three simulations settings comparable to MDE, performing particularly well in the experimental sample.  The remaining methods struggle in at least one of the comparisons.



% latex table generated in R 3.3.1 by xtable 1.8-2 package
% Mon Jan 16 11:58:45 2017
% latex table generated in R 3.3.1 by xtable 1.8-2 package
% Mon Jan 16 19:41:59 2017
\begin{table}[tb!]
\centering\spacingset{1}
\begin{tabular}{r|rr:rr}
  \hline
Outcome&  \multicolumn{2}{c}{Observational Untreated}&  \multicolumn{2}{:c}{Experimental Untreated}\\
 & Bias & RMSE & Bias &  RMSE \\
  \hline
OLS, in-sample & 0.00 & 10173.27 & 0.00 & 5433.39 \\
\hdashline
 MDE & -15054.61 & 21267.99 & 428.78 & 7040.68 \\
  MDE, bagged& -255.94 & 10804.69 & 660.84 & 6865.75 \\
  Horseshoe& -16104.93 & 22371.13 & 14733.41 & 15798.17 \\
  BART & -11049.12 & 17046.63 & 1125.95 & 6896.35 \\
  POLYMARS & -13035.59 & 18957.33 & 996.27 & 6646.84 \\
  SuperLearner & -13407.01 & 19474.82 & 2057.12 & 6684.54 \\
   LASSO & -12942.43 & 19186.10 & 1689.17 & 6751.50 \\
   \hline
   \hline
\end{tabular}\caption{{\bf Prediction Exercise, LaLonde Data.}  Columns contain the results from predicting the outcome on held-out subsets of the LaLonde data (top row).   We fit a model to the experimental data and then use this model to predict outcomes in the observational data (columns 2-3).  Next, we fit a model to the experimental treated and observational untreated, then use this model to predict the held-out experimental group (columns 4-5).  In each case, we include the result from least-squares fit to the held-out sample as a baseline.  We see that the bagged MDE estimate performs well in both settings, in terms of both bias and RMSE.  In the second setting, POLYMARS achieves the lowest RMSE, but at the cost of  a  large bias.}\label{tab:lalonde2}
\end{table}



Throughout this study, we have focused not just on sample-average estimates but individual-level effect estimation.  We next compare several methods on their ability to predict a held-out subsample. Results are presented in Table \ref{tab:lalonde2}. In the first comparison (columns 2-3), we fit a model to the experimental data and then use this model to predict outcomes in the observational data.  As an experimental sample may not resemble a target population of interest, this prediction exercise indicates the extent to which methods can generalize from an experiment to observational data.  In the second comparison (columns 4-5), we fit a model to the experimental treated and observational untreated, then use this model to predict the held-out experimental group.  In each case, we include the result from least-squares fit directly to the held-out sample as a baseline.

We see again that bagged MDE performs well in both exercises, particularly the first, where it achieves a predictive RMSE nearly that of OLS fit to the held-out data.  In the first simulation, MDE performs less well, achieving a higher RMSE than all but the horseshoe.  In the second analysis, predicting the experimental untreated, bagged MDE achieves a RMSE and bias lower than BART.  MDE achieves the lowest bias, but with a somewhat higher RMSE than the remaining methods except for the horseshoe.  POLYMARS and the SuperLearner both achieve a lower RMSE than MDE and bagged MDE in the second analysis, but at the cost of a higher variance.




\subsection{Application: Instrumental Variables}

Next we examine an empirical application using instrumental variable methods. Here we focus on an application with a continuous instrument in order to highlight how the proposed methodology naturally incorporates continuous as well as non-continuous exogenous variables.

\citet{larreguy2017} study the long term political effects of increased education. To do this they utilize variation in the intensity of a Nigerian government reform, the Universal Primary Education reform of 1976. The authors leverage Afrobarometer survey data to explore a variety of political variables, such as interest in the news and knowledge of politics. Concerned about endogeneity problems present in regressing these political variables onto measures of education, they use an instrumental variables strategy akin to earlier work in development economics on the effects of education \citep{duflo2001}. In particular they exploit temporal (the program started in 1976 and so impacted particular cohorts of citizens) and spatial (regional level differences between actual and potential enrollment, see also \citet{bleakley2010malaria}) variation. They use linear two stage least squares with an interaction between the intensity of the program and whether an individual would have been eligible to benefit from the program as an instrument for education levels. The authors include a range of control variables and fixed effects for cohort, region, and Afrobarometer survey wave. They find strong and robust impacts of education on long term political variables. Using their replication data we implement our proposed method to reevaluate their results.

Part of our interest is in allowing for the instrument to have non-linear effects on the endogenous variable. Hence before displaying results from our analysis, we present evidence of such a non-linear relationship by simply fitting a generalized additive model to the data, allowing for local smoothing on the instrument and including the full set of controls. On the x-axis we present the values of the instrument along with a histogram of its marginal distribution. On the y-axis we present the fitted values of the endogenous variable.\footnote{We re-scaled the endogenous variable, education, which in the original analysis ran from 1 to 5.} Here and below we present this distribution as a 2-dimensional heatmap in order to convey density of observations over the space, and fit a smooth trend line to convey the basic pattern. Generally speaking we see higher effects of the instrument at lower values of the instrument, and lower effects at higher values of the instrument. This gives some evidence that it might be desirable in calculating causal effects to incorporate greater functional form flexibility into the analysis.

\begin{figure}[t!]
%\vspace{-.9in}
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{GAM-Instrument.pdf}
    \vspace{-.25in}
  \end{center}\spacingset{1}
  \caption{{\bf Nonlinear bi-variate relationship between instrument and fitted values of the endogenous variable. Estimated using generalized additive model and not the proposed method.}      } \label{fig:GAM1}
\end{figure}

We fit our model using a full set of splines and interactions described in Section~\ref{sec:estimation}.\footnote{The original model used a substantial number of fixed effects, in part to establish a difference in difference identification strategy with an instrumental variable. This poses no problem for our proposed method. However for computational purposes we made some minor modifications that we discuss in Appendix~\ref{app:simdata}.} We present two main sets of results.

First, in Figure~\ref{fig:firststage} we present a plot analogous to Figure~\ref{fig:GAM1} which plots for each observation the relationship between the instrument and the covariance between the instrument and and (endogenous) treatment variable. This is the ``first stage'' result. We see a strong positive relationship until the upper end of the instrument distribution, where it begins to get closer to 0. A behavioral interpretation of this pattern is that the returns on education to increasing intensity in the program were lower in areas with greater intensity, perhaps because it would be harder to increase their education levels even higher.\footnote{We also investigated the first stage fit of our model compared to least squares via cross-validation. We found nearly identical performance, which is impressive for our proposed method given that the sample size is quite large relative to the number of parameters fit by the least squares model.}

\begin{figure}[t!]
%\vspace{-.9in}
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{InstrumentVersusFirstStage.pdf}
    \vspace{-.25in}
  \end{center}\spacingset{1}
  \caption{{\bf Relationship between instrument and first stage covariance (covariance between instrument and treatment). Estimated using the proposed method.}      } \label{fig:firststage}
\end{figure}

Figure~\ref{fig:secstage} plots the ``second stage'': the treatment variable versus the covariance between the treatment and outcome variable. A piecewise line plots the means across the values of the treatment variable. At lower lower levels of the treatment variable the average covariance was lower than at higher levels of the treatment variable. Combined with the results in Figure~\ref{fig:firststage} this helps to give us an expectation of what the LICE will look like over the sample.
%iv1$scaled.2<-(iv1$mod2$cf.fit-iv1$mod2$fit)/0.001


\begin{figure}[t!]
%\vspace{-.9in}
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{TreatmentVsScaled2ndStage.pdf}
    \vspace{-.25in}
  \end{center}\spacingset{1}
  \caption{{\bf Treatment variable versus covariance between treatment and outcome.}} \label{fig:secstage}
\end{figure}


Remembering that the LICE is a ratio estimate, the next step then is to directly plot the numerator of our IV estimate (the covariance between the instrument and treatment) against the denominator (the covariance between the treatment and outcome).  Figure~\ref{fig:firstvsec} plots the results. Several patterns are interesting. First, the majority of observations are in the upper right quadrant, with positive covariances. Second, there is a non-linear relationship. As the first stage covariance increases (that is, the strength of the instrument), there are declining returns to having an effect on the second stage relationship. Substantively this implies a positive relationship overall between education and interest in news, but one that is diminishing at higher levels once endogeneity concerns have been addressed. Third, there are a small number of observations in the bottom left quadrant. For these observations the LICE will still be positive. On average these observations are likely to be of lower education (see Figure~\ref{fig:secstage}). Behaviorally these individuals were negatively encouraged but also saw a decrease in interest in news. Finally, there are no observations in the top left quadrant, and few observations in the bottom right. Individuals in the bottom right were positively encouraged by the reforms but saw a decrease in news interest; for these individuals we would expect a negative LICE

\begin{figure}[t!]
%\vspace{-.9in}
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{FirstStageVs2ndStage.pdf}
    \vspace{-.25in}
  \end{center}\spacingset{1}
  \caption{{\bf First stage versus second stage.}} \label{fig:firstvsec}
\end{figure}




Finally in Figure~\ref{fig:ICE} we plot the distribution of local individual causal effects (LICE) against the treatment. The dashed horizontal lines represents the average LICE, which we estimate to be .35 (95\% confidence interval: .20, .49).\footnote{We bootstrapped the confidence intervals using 100 bootstrap runs. Using TSLS the original analysis returned a point estimate of .62. \citet{larreguy2017} cluster standard errors at the state level. However, re-analysis of their specification shows this made no difference. In fact, the confidence intervals were slightly wider without clustering. We did not cluster standard errors, though a block bootstrap approach would accomplish this goal.} While there was non-linearity in the first stage estimate, we observe more stability over the sample in the LICE. This directly follows from the results in Figure~\ref{fig:firstvsec}.


\begin{figure}[t!]
%\vspace{-.9in}
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{TreatmentVersusICE.pdf}
    \vspace{-.25in}
  \end{center}\spacingset{1}
  \caption{{\bf Local Individual Causal Effect versus treatment variable.}      } \label{fig:ICE}
\end{figure}









\section{Conclusion}
\label{sec:conclusion}


This paper has focused on using a high-dimensional, flexible regression model directly target observation-level causal effect estimates by considering an extremely rich covariate \emph{and} functional space. We utilize recent advances in machine learning to demonstrate that the proposed approach performs extremely well against other cutting edge methods. Finally, we show how our estimation strategy extends beyond the estimation of simple treatment effects and can also be used in instrumental variable or other structural models.

%This paper has focused on using a high-dimensional, flexible regression model to estimate observation level causal effects, rather than targeting a low-dimensional structural parameter. While a core building block of the causal inference literature, highlighted by the use of potential outcomes, most statistical tools used by applied researchers target a structural parameter. We instead directly target individual level causal effect estimates by considering an extremely rich covariate \emph{and} functional space. We utilize recent advances in machine learning to demonstrate that the proposed approach performs extremely well against other cutting edge methods. Finally, we should how our conceptualization easily extends beyond the estimation of simple treatment effects, but can also be used in encouragement designs (instrumental variables) and the analysis of causal mechanisms.


%In many ways, this work is motivated by estimation, but the close focus on estimation brought several causal insights to the fore.  We build off  seminal work on  multivariate spline models \citep{Stone1997,friedman1991multivariate} but combine it with more recent machine learning tools \citep{Fan:2008,RatkovicTingley17}.  The reliability and flexibility of the model gives us two advantages.  First, we can sidestep inverse weights, matching, subclassifying, etc., which are rely on a treatment model that is rarely of interest.  The focus on the outcome model allows us to devote our energy to getting the outcome model correct.  We show above that this is a fruitful line of inquiry.  Second, the estimation strategy allowed us to formulate and estimate individual-level counterfactual estimates.  Doing so allows us to combine the observation-level potential outcomes approach of Rubin with Pearl's nonparametric SEM framework in a natural manner.

We see this work as providing an integrative thread across a number of seminal contributions. Our approach starts with Rubin's observation that causal inference is a missing data problem: were the counterfactual outcomes known, causal inference would amount to simple descriptive statistics.   We work towards this goal head-on, avoiding workarounds like inverse weights and propensity scores, which add another layer of modeling uncertainty and are rarely subjects of direct interest.  An early stumbling block in causal inference was acknowledged by Robins, in that causal inference is only as persuasive as the underlying model.  Rather than develop methods that are reasonable even in the face of model misspecification, we focus our attention on getting the conditional mean correct.  We build off  seminal work on  multivariate spline models \citep{Stone1997,friedman1991multivariate} but expand the range of basis functions utilized and combine it with recent machine learning tools \citep{Fan:2008,RatkovicTingley17} in order to focus on this goal.\footnote{Of note, as shown in Section~\ref{sec:simulation}, our expanded basis function approach also enabled the LASSO to beat a number of cutting edge machine learning methods.}  Though we started the project in terms of the potential outcomes approach, we found that our estimation strategy combined thinking in terms of observation-level counterfactuals as well as structural ``do'' manipulations, as developed by Pearl.  The utility of our approach becomes more clear in the case of the instrumental variable analysis, where we integrate the nonparametric structural equation model (SEM) approach and the potential outcome approach of \citet{angrist1996identification}.



The are a number of areas for future work. One thing we have left out of the current paper is the set of advantages of using the Bayesian LASSOPlus that we previously developed \citep{RatkovicTingley17}.\footnote{These advantages include straightforward ways to incorporate binary or truncated outcome variables, random effects, uncertainty estimates with desirable coverage properties, and not relying on arbitrary tuning parameter selection.} We are actively extending our approach to cases with multiple separate treatment variables, instruments \citep[e.g.,][]{chernozhukovinstrumental,belloni2012sparse}, or mediators. Our current approach considers the impact of perturbing a single instrument, but extending it to multiple instruments and treatments will involve moving from partial to Frechet derivatives.  Future work could also consider how to extend the model to the longitudinal setting or how to incorporate spline bases with discontinuities so as to capture ``jumps'' in the data, and how to estimate the effect in a regression discontinuity design.

\clearpage
\bibliography{arxiv_bibliography}
\clearpage

\appendix

\clearpage
\spacingset{1.5}


\Large \textbf{Appendix}\\\bigskip
\normalsize
The outline of the appendix is as follows.\\\medskip
Section~\ref{app:lassoplus} details properties of the LASSOPlus estimator. Interested readers can consult \citet{RatkovicTingley17} for additional details.  Section~\ref{app:bootstrap} proves validity of our bootstrap estimate.  Section~\ref{app:ivsims} conducts a simulation study on the IV compliance threshold and proposes a diagnostic for instrumental variables akin to the first-stage $F$-statistic. Section~\ref{app:proofs} gives the proofs of our Oracle results. Section~\ref{app:simdata} discusses our processing of the \citet{larreguy2017} data as well as plots variable effect estimates.






\section{LASSOplus Properties}\label{app:lassoplus}


Our joint conditional density on $c, w$ produces a problem similar to the adaptive LASSO of  \citet{zouadaptive}:
\begin{gather}
- \log\l(\Pr\l(c, \{w_k\}_{k=1}^p|\lambda, \sigma^2, \gamma\r)\r) =\frac{1}{\sigma^2}\l\{ \frac{1}{2} \sum_{i=1}^n (Y_i-R(T_i, X_i)^\top c)^2 +\lambda \sigma\sum_{k=1}^p w_k |c_k|\r\}+\sum_{k=1}^p w_k^{\gamma}.
\end{gather}
where the weights and global sparsity parameter enter the log-posterior as $w_k^\gamma$ rather than as separate terms, as in  \citet{Leng2014,  Alhamzawi2012, Griffin2012, Griffin2010, Kang2009}.



\paragraph{Estimation} Following PC, we reintroduce conjugacy in the mean parameters through augmentation:
\begin{gather}
c_k \sim DE(w_k \lambda/\sigma)  \Rightarrow c_k|\tau_k^2, \sigma^2 \sim \mathcal N(0, \tau_k^2\sigma^2); \;\;\tau_k^2\sim \exp(\lambda^2w_k^2/2).
\end{gather}
 Both the MCMC and EM estimation details are standard, see \citet{RatkovicTingley17}.

\paragraph{The tuning parameter.}  We first focus on the conditional posterior density of the tuning parameter, $\lambda$.  Rescaling by $n$ reveals that we are recovering the MAP estimate
\begin{gather}
\widehat c\;|\cdot = \argmin_c \frac{1}{n} \sum_{i=1}^n \frac{(Y_i-R(T_i, X_i)^\top c)^2}{2\sigma^2} +\frac{\lambda}{ n \sigma}\sum_{k=1}^p w_k |c_k|.\label{eq:map1}
\end{gather}
The Oracle rate is achieved when $\lambda/n$ grows as $\sqrt{\log(p)/n}$, i.e. when $\lambda$ grows as $\sqrt{n \log(p)}$.  Our conditional posterior density of $\lambda^2$ is
\begin{align}
\lambda^2|\cdot \sim \Gamma\l(n\times\l(\log(n) + 2 \log (p) \r) , \sum_{k=1}^p \tau^2_k/2+\rho\r)\label{eq:lambdacons}\\
\end{align}
which several desirable properties.  First, each mean parameter $c_k$ is given its own shrinkage parameter, $\frac{\lambda w_k}{n\sigma}$, allowing for differential regularization of large and small effects.  Second, for $p \sim n^\alpha, \alpha>0$, then $\widehat \lambda =O\l(\sqrt {n \log(p)} \r)$, as desired.  Second, the tuning parameter has the structure  given in \citet[][Corollary 6.2]{Buhlmann2013} so that it provides a consistent estimate of the conditional mean.  The term $\log(n)$ controls the probability with which the Oracle bound holds, and it goes to zero in $n$.


\paragraph{The weights.}
To see the impact of the weights, consider the posterior conditional density and mean of $w_k$,
\begin{align}
\Pr(w_k|\cdot) =\frac{ e^{-w_k^\gamma - \frac{\lambda w_k}{\sigma}|c_k|}}{ \int e^{-w_k^\gamma - \frac{\lambda w_k}{\sigma}|c_k|} dw_k}; \quad \E(w_k|\cdot) = \int_{w_k=0}^\infty w_k \Pr(w_k|\cdot) dw_k
\end{align}
with $\widehat w_k$ simply the conditional mean of $w_k$ under this density.

The derivative of $\widehat w_k$ with respect to $|c_k|$ evaluated at the EM estimate gives
\begin{align}
\frac{d \widehat w_k}{d|c_k|}=-\widehat \lambda \sqrt {\widehat{\frac{1}{\sigma^2} }}\Var(w_k|\cdot)<0
\end{align}
showing that the weights are inversely related to the magnitude of the estimates, $|c_k|$.   This inverse relationship between the weights and the estimate gives us the same properties as the adaptive LASSO of \citet{zouadaptive}.


Next, we consider the two limiting cases $|\widehat c_k| =0$ and $|\widehat c_k|\rightarrow \infty$.  Consider first $|\widehat c_k|\rightarrow \infty$.  In this case, the exponent in the conditional kernel of $w_k$ is dominated by the term $-\frac{\lambda w_k}{\sigma}|c_k|$ and approaches an exponential density with mean $\widehat w_k \rightarrow \widehat \sigma/\widehat \lambda$.  Therefore, as $|c_k|$ grows, the weighted LASSO penalty approaches $1/n$, a negligible term.

When $|\widehat c_k|=0$, the conditional posterior density follows a generalized Gamma density with kernel $\exp\{-w_k^{\widehat \gamma}\}$ and has mean $\Gamma(2/\widehat \gamma)/\Gamma(1/\widehat \gamma)$ which we denote $\overline \gamma$.\begin{footnote}{ $\Gamma(\cdot)$ refers to the gamma function (not density)}\end{footnote}  Therefore, for the zeroed out parameters, the weighted LASSO penalty for $|c_k|$ approaches $\widehat \lambda \overline \gamma/n$.  Since $\widehat \lambda = O(\sqrt{n \log(n)})$, the penalty is of order $\log(n)/\sqrt n$, which is not negligible.  The value of this parameter is controlled by $\widehat \gamma$, to which we turn next.

\paragraph{The global sparsity parameter.}


The global sparsity parameter $\gamma$ serves to pool information across the mean parameters \citep[see e.g.][for a similar insight]{bhattacharya2015dirichlet}.  To illustrate its workings, we consider the two limiting cases, $\gamma \in \{0, \infty\}$.

First, as $\gamma$ approaches $0$, our prior approaches the spike-and-slab prior for a mean parameter, resulting in no shrinkage but with a point mass at zero.  Taking $\gamma \rightarrow 0 \Rightarrow \Pr(w_k|\cdot) \rightarrow Exp(\lambda|c_k|/\sigma)$.  This gives us $\widehat w_k =  \l(\widehat{\frac{1}{\sigma^2} }\r)^{-1/2}/(\widehat \lambda |\widehat c_k|)$.  Plugging into the prior on $c_k$ gives $\Pr(c_k)\rightarrow DE\l(\widehat \lambda \widehat w_k   \l(\widehat{\frac{1}{\sigma^2} }\r)^{1/2} \r)\rightarrow DE\l(\widehat \lambda \l(\widehat{\frac{1}{\sigma^2} }\r)^{-1/2}/(\widehat \lambda |\widehat c_k|)  \l(\widehat{\frac{1}{\sigma^2} }\r)^{1/2} \r) \propto 1$, the flat Jeffreys prior when $\widehat c_k \ne 0$.  When $\widehat c_k=0$, the prior has infinite density, to due the normalizing constant of order $1/|\widehat c_k|$, giving a spike-and-slab prior with support over the real line.

Taking $\gamma \rightarrow \infty \Rightarrow \Pr(w_k|\cdot) \rightarrow U(0,1)$, a uniform on $[0,1]$, since any weight greater than $1$ has mass proportional to $\exp\l\{w^{-\gamma}\r\}=0$.  This gives us $\widehat w_k = 1/2$.  Plugging into the prior on $c_k$ gives $\Pr(c_k)\rightarrow DE\l(\frac{1}{2}\widehat \lambda   \l(\widehat{\frac{1}{\sigma^2} }\r)^{1/2} \r)$, which is the PC prior with $1/2$ the rate parameter.

\section{Proof of Validity of Bootstrap Estimate}\label{app:bootstrap}

Denote as $T_{B_N}$ the complete set of all possible bootstrap estimates.
We need to show that the bootstrap estimate $T_{B_N}$ converges in distribution to $\widetilde T_i$.  We exploit the independence between the bootstrapped fitted values, $\widehat T_{B_N}$ and permuted bootstrapped errors $\widehat \epsilon_{B^\prime_N}$, as well as Slutsky's Theorem, to show that a bootstrapped approximation to the convolution formula converges to its population analog. We assume that the distribution of $\widetilde T_i|X_i$  is Glivenko-Cantelli with bounded density a.e. and that  $\epsilon_i \indep X_i$.  Then, with $F_{B_N}$ the empirical bootstrap distribution over replicates
\begin{align}
F_{B_N}(z) &= \Pr(T_{ib} \le z)=\Pr(\widehat T_{ib} + \widehat \epsilon_{ib^\prime} \le z) \\
&=\sum_{i=1}^{B_N} \sum_{x\in \widehat \epsilon_{B^\prime_N}}  \frac{1}{B_N} \mathbf 1(x=\widehat \epsilon_{ib}) \mathbf 1(\widehat T_{ib} + x \le z)   \convp \int f_{\epsilon_i}(x) F_{\widetilde T_i} (z-x)dx\\
&=F(\widetilde T\le z)
\end{align}
by Glivenko-Cantelli and Slutsky.  In practice $T_{B_N}$ is too large to compute, so we use $T_{B}$ with large $B$.

\section{Instrumental Variables Diagnostics and Simulations}
\label{app:ivsims}


\subsection{Simulations}


We next present a short set of simulation results to examine the performance of our method applied to the case of an encouragement design (instrumental variable). For simplicity we compare our proposed method to two-stage least squares estimation.\footnote{The conclusion discusses connections with \citet{chernozhukovinstrumental,belloni2012sparse} that investigate the case of variable selection, rather than variable and functional selection, in the estimation of a structural parameter rather than individual effects. Methods used in Section~\ref{sec:simulation} have not been extended to the instrumental variables context except for ensemble tree methods by \citet{athey2016solving}.} We examine first stage and second stage performance, as well as how well the proposed method helps to unpack individual level compliance estimates




\paragraph{Simulation Environments}  We consider five different simulation settings.   We fix the second stage (outcome) model but vary the first stage (treatment) model across five settings.  The fist is a null model, in which the instrument does not impact the treatment and therefore no observation complies with the instrument.  In the second, the instrument has an additive linear effect on the treatment.  The third scenarios contains a mixture of compliers, non-compliers, and defiers.  The fourth scenario is the non-linear model from the earlier simulation setting, and the fifth contains a discontinuity in the instrument. We note that, in this last setting, the true function is not in the model space of our spline bases.


Specifically, we use the exact same settings as in the previous simulations to generate the treatment from the instrument.  We add an additional null model,
\begin{align}
\textrm{Null }&\;T_i=\sum_{k=5}^8 X_{ik} \beta_k+\epsilon_i^T,
\end{align}
resulting in five first stage models. The instrument is also generated some confounding from covariates
\begin{align}
Z_i =-3+(X_{i1}+X_{i4})^2+\epsilon_i^T; \;\;\epsilon_i^Z\iid \mathcal N\l(0, 4\r).
\end{align}
Across settings, the outcome is generated as
\begin{align}
Y_i = |T_i-2|+\sum_{k\in\{5,6,7,8\}} X_{ik} \beta_k +X_{i1} \times X_{i2} +\epsilon_i^Y.
\end{align}
where
\begin{align}
[\epsilon_i^T, \epsilon_i^Y]^\top \iid N\l([0,0]^\top, C \begin{bmatrix}1, .9\\.9,1 \end{bmatrix}\r)
\end{align}
where $C$ is selected such that the second stage has a signal to noise ratio of 1. We generated the instrument similar to above, with nonlinear confounding with two pre-treatment covariats.  The error structure was induced so endogeneity bias would be severe.




\paragraph{Results}



In each setting, we compare the performance of MDE to TSLS.  We consider three sets of performance measures. The first two are RMSE performance for the first and second stages.  The second stage RMSE is only measured off the observations with an in-truth identified causal effect.  Third, we consider the ability of MDE to differentiate compliers from defiers and for our threshold to identify non-compliers.


We focus on two ways to evaluate our results. First, we examine the RMSE on the first and second stage. The RMSE in the second stage is particularly important because this evaluates the extent to which we are capturing the LICE. Results for the RMSE in the first stage are still helpful to show, though they are analogous to the insights provided in Section~\ref{sec:simulation}.


Figure~\ref{fig:iv_ratio} plots the ratio of the proposed method's RMSE to the RMSE of 2SLS. Values less than one indicate better performance for the proposed method.\footnote{Figures~\ref{fig:iv_first} and ~\ref{fig:iv_second} in Appendix~\ref{app:simdata} plot the actual RMSE's for both methods.} In all settings except for the linear model, and for both the first and the second stage, the method of direct estimation returns a superior RMSE. As the sample size increases in the linear model, MDE quickly catches up.

\begin{figure}[t!]
%\vspace{-.9in}=
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{IVSims-RMSERatio.pdf}
 %   \vspace{-2.2in}
  \end{center}\spacingset{1}
  \caption{{\bf Ratio of RMSE estimates for proposed method versus 2SLS for first stage and second stage estimation across simulation environments.}      } \label{fig:iv_ratio}
\end{figure}



\paragraph{Compliance Estimates}

Figure~\ref{fig:iv_compliance} examines how well we estimate the correct sign for first stage. In our simulations we know if someone was positively encouraged by the instrument, negatively encouraged, and not encouraged (zero).\footnote{Figure~\ref{fig:iv_compliance_true} in Appendix~\ref{app:simdata} gives the true proportions in the simulated data. } We then recorded the model's individual level estimates and calculate the percentage of observations correctly classified into each bin. We make several observations on the figure.  First, where there are in-truth no observations encouraged, as in the first row, the proposed threshold does indeed return this value. In the second row, the first stage estimate is an additive linear model in which all observations comply positively, and again, the threshold separates compliers from others even at a modest sample size. In the interactive setting the performance is increasing in sample size for identifying both positive and negative compliers. In the non-linear setting our performance for positive observations increases in sample but the proposed method performs poorly for the other types. Finally in the discontinuous setting positive and negative observations are increasingly better identified. Zero types are estimated with high accuracy, though there is a slight degradation as the sample increases. %Negative observations are recovered at a low rate but this performance increases slowly in the sample size.\footnote{In ongoing work we are working to improve the performance  of the threshold  by studying how misclassification is happening and how this relates to our analytically derived threshold.}


We also recorded compliance estimates for TSLS.  TSLS can only return one compliance estimate for each observation, and we recorded all observations as not complying if the first-stage $F$ statistic on the instrument was less than $10$.  Results are in Figure~\ref{fig:iv_first} of Appendix~\ref{app:simdata}. When TSLS did recover compliance percentages correctly, it was only for positively encouraged units. Given the large number of other types in the simulations this performance is not desirable.

\begin{figure}[t!]
%\vspace{-.9in}=
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{IVSims-CompliancePercentages-Stage1.pdf}
 %   \vspace{-2.2in}
  \end{center}\spacingset{1}
  \caption{{\bf First stage accuracy for proposed method.} Blue line plots the percentage of in fact positive LICE that we correctly capture. The red line plots the percentage of in fact negative LICE that we correctly capture. The green line plots the percentage of in fact non-identified units that we correctly capture. In setting 1 no one is identified. In our simulations there are no observations with in truth 0 LICE and so the purple line for this case is not present.} \label{fig:iv_compliance}
\end{figure}





\subsubsection{True Compliance Rates by Simulation Setting}



Figure~\ref{fig:iv_compliance_true} presents the true percentage for the first stage compliance status. This shows how in setting 1 no observations were encouraged. Settings 2-5 mix the ratio in different ways across positive, negative, and in truth 0.


\begin{figure}[t!]
%\vspace{-.9in}=
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{IVSims-CompliancePercentages-1-True.pdf}
 %   \vspace{-2.2in}
  \end{center}\spacingset{1}
  \caption{{\bf Percentage of observations in sample by compliance status.}      } \label{fig:iv_compliance_true}
\end{figure}





Figure~\ref{fig:iv_compliance_true} presents the true percentage each sign for the instrumented causal effect of observations. This shows how in setting 1 all observations do not have an identified sign in the second stage. Settings 2-5 mix the ratio in different ways across positive, negative, in truth 0, and unidentified.


\begin{figure}[t!]
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{IVSims-CompliancePercentages-True.pdf}
  \end{center}\spacingset{1}
  \caption{{\bf Percentage of observations in sample by category of sign on instrumented causal effect.}      } \label{fig:iv_compliance_true2}
\end{figure}


\subsubsection{RMSE}

Figure~\ref{fig:iv_first} presents RMSE estimates for the first stage model. Consistent with the previous simulation results, the proposed method performs well across a range of settings. The standard least squares model, which is misspecified, generally does not perform as well, especially in the third simulation environment.


\begin{figure}[t!]
%\vspace{-.9in}=
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{IVSims-FirstStage.pdf}
 %   \vspace{-2.2in}
  \end{center}\spacingset{1}
  \caption{{\bf RMSE estimates for proposed method and 2SLS for first stage estimation across simulation environments.}      } \label{fig:iv_first}
\end{figure}

Figure~\ref{fig:iv_secon} presents the second stage results. We again see improved performance by the proposed method. RMSE's are lower than 2SLS estimates. In the simulation context we consider there are mixtures of individuals compliance types (discussed further below). Not capturing this heterogeneity means our estimate of the LICE is more accurate than the 2SLS.



%Across all but the first simulation environment the proposed method has a higher RMSE at low sample sizes. While this might seem at odds with the poorer performance of the standard 2SLS first stage specification, the good second stage performance is because X.

%https://pdfs.semanticscholar.org/ab95/86b51f8685abc347ceb4619a97a4b5586d9b.pdf

\begin{figure}[t!]
%\vspace{-.9in}=
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{IVSims-SecondStage.pdf}
 %   \vspace{-2.2in}
  \end{center}\spacingset{1}
  \caption{{\bf RMSE estimates for proposed method and 2SLS for second stage estimation across simulation environments.}      } \label{fig:iv_secon}
\end{figure}



\subsubsection{TSLS Peformance}

Next we give compliance estimates for TSLS in the first stage.  TSLS can only return one compliance estimate for each observation, and we recorded all observations as not complying if the first-stage $F$ statistic on the instrument was less than $10$.  Results are in Figure~\ref{fig:iv_first}. When TSLS did recover compliance percentages correctly, it was only for those with a positive instrumented causal effect. Essentially it records everything as falling in this category. Given the large number of other types in the simulations this performance is not desirable.

\begin{figure}[t!]
%\vspace{-.9in}=
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{IVSims-CompliancePercentages-Stage1-IV.pdf}
 %   \vspace{-2.2in}
  \end{center}\spacingset{1}
  \caption{{\bf First stage accuracy for TSLS. Blue line plots the percentage of in fact positive LICE that TSLS correctly captures. The red line plots the percentage of in fact positively encouraged that TSLS correctly captures. The green line plots the percentage of negatively encouraged units that TSLS correctly captures and the blue line captures the percentage of non-encouraged observations correctly identified.}      } \label{fig:iv_first}
\end{figure}


We also recorded compliance estimates for TSLS in the second stage.  TSLS can only return one compliance estimate for each observation, and we recorded all observations as not complying if the first-stage $F$ statistic on the instrument was less than $10$.  Results are in Figure~\ref{fig:iv_secon}. When TSLS did recover compliance percentages correctly, it was only for those with a positive instrumented causal effect. Given the large number of other types in the simulations this performance is not desirable.

\begin{figure}[t!]
%\vspace{-.9in}=
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{IVSims-CompliancePercentages-2SLS.pdf}
 %   \vspace{-2.2in}
  \end{center}\spacingset{1}
  \caption{{\bf Second stage accuracy for TSLS. Blue line plots the percentage of in fact positive LICE that TSLS correctly captures. The red line plots the percentage of in fact negative LICE that TSLS correctly captures. The green line plots the percentage of in fact non-identified units that TSLS correctly captures.}      } \label{fig:iv_second}
\end{figure}



\subsection{Diagnostics}

Instrumental variables estimation is reliable to the extent that the encouragement of the instrument has a strong effect on the treatment \citep{bound1995problems,staiger1994instrumental}.  The standard measure of this effect in the TSLS case is the $F$-statistic.\footnote{Though see \citet{kleibergen2002pivotal} for one alternative proposal.} Many researcher use a threshold, like 10, to decide whether their instrument is strong enough.

Using the notation in the paper, the first-stage $F$-statistic for our proposed method is:
\begin{align}
\widehat F = \frac{\frac{1}{\widehat{df}^\nabla}\sum_{i=1}^n \l\{R_T^{ZX}(Z_i, X_i) ^\top \widehat c_T^{ZX} - \widehat \mu_T \r\}^2}{\frac{1}{n-\widehat df^\nabla}\sum_{i=1}^n \l\{T_i - R_T(Z_i, X_i) ^\top \widehat c_T  - \widehat \mu_T\r\}^2}
\end{align}
where $R_T^{ZX}(Z_i, X_i)$ is the the subvector of $R_T(Z_i, X_i)$ that fluctuates with $Z_i$ and $c_T^{ZX}$ the corresponding parameters.




\section{Oracle Proofs}\label{app:proofs}

We condition on $\widehat W=diag(\widehat w_k)$ throughout the proof and note that $\widehat W=I_p$ reduces the proof to that of the standard LASSO.  Throughout, to ease notation, we write $R=R(T_i, Z_i)$, etc. dropping dependence on $T_i, Z_i$ when it is clear from context.  We also assume that $Y$ and $R_k$ have sample mean zero, so we do not have to worry about an intercept.


Start with the LASSO problem
\begin{align}
\widehat c = \argmin_{\widetilde c} ||Y-R \widetilde c||_2^2 + \lambda ||\widehat W \widetilde c||_1^1.
\end{align}
As we are minimizing over the sample, we know the estimator $\widehat c$ satisfies
\begin{align}
 ||Y-R \widehat c||_2^2 + \lambda ||\widehat W  \widehat c||_1^1 \le ||Y-R c||_2^2 + \lambda || W c||_1^1.
\end{align}
After some manipulation \citep[e.g.,][]{RatkovicTingley17, Buhlmann2013}, this excess risk generates the inequality
\begin{align}
\frac{1}{n}\l\{||R \widehat \delta||_2^2+ \lambda||\widehat W \widehat \delta||_1^1\r\}\le
C\frac{\widehat \lambda^2 \widehat \sigma^2 \overline \gamma^2 |S|}{n^2 \phi_0^2}+ C_\infty \frac{||R^o_{\infty/n}c_{\infty/n}||^2_\infty}{n} +C_\perp\frac{  ||R^o_{\perp}c_{\perp}||_\infty^2 }{n}
\end{align}
which is Equation \ref{eq:oracle1}.

\paragraph{First order conditions on the LASSO and the derivative.}  The first-order conditions for the LASSO problem above are
\begin{align}
2 R_k^\top  \widehat \epsilon=s_k \widehat \lambda \widehat w_k \label{eq:opt1}
\end{align}
with $s_k \in[-1,1]$ and $s_k =1$ or $-1$ iff $\widehat c_k \ne 0$.  We consider the subset of the matrix corresponding with these non-zero estimates
\begin{align}
R^{\widehat S}=\{R_k: \widehat c_k \ne 0\}
\end{align}
and the submatrices $R^{\widehat S, X}$ and $R^{\widehat S, XT}$.


Consider the gradient derivative of equality \ref{eq:opt1} wrt the treatment at each observation and noting, by the identification assumptions (1)-(4), $\frac{\partial}{\partial T_i} \widehat s_k \widehat w_k \widehat \lambda=0$, which gives
\begin{align}
0&=\sum_{i=1}^N \frac{\partial}{\partial T_i} \sum_{i=1}^N  \l(R_{ik}^{\widehat S} \widehat \epsilon_i\r) \Rightarrow \\
0 &=\sum_{i=1}^N \frac{\partial}{\partial T_i} R_{ik}^{\widehat S} \widehat \epsilon_i\\
&= \sum_{i=1}^N \frac{\partial R_{ik}^{\widehat S} }{\partial T_i} \widehat \epsilon_i+ R_{ik}^{\widehat S} \frac{\partial \widehat \epsilon_i}{\partial T_i}
\end{align}
This gives us
\begin{align}
\bigg|\sum_{i=1}^N \frac{\partial R_{ik}^{\widehat S} }{\partial T_i} \widehat \epsilon_i\bigg| = \bigg|\sum_{i=1}^n  R_{ik}^{\widehat S} \frac{\partial \widehat \epsilon_i}{\partial T_i}\bigg|
\end{align}
which is a solution to the following LASSO problem
\begin{align}
\widehat c^{S,T} = \argmin_{c^{S,T}} ||\widetilde{\Delta Y}-\Delta R^{\widehat S,T} c^{\widehat S,T}||_2^2 + ||\widetilde W^{\widehat S, T} c^{\widehat S,T}||_1^1
\end{align}
where $\Delta R^{\widehat S,T}$ is the elementwise partial derivative of $R^{\widehat S,T}$ wrt $T_i$ and
 $\widetilde{\Delta Y}$ is a pseudo-response constructed from the estimate $\widehat{\Delta Y}=\l[\widehat{\frac{\partial Y_i}{\partial T_i}}\r]$ as
\begin{align}
\widetilde{\Delta Y} &=\Delta R^{\widehat S,T}\widehat c^{\widehat S,T}+\widehat \epsilon.
\end{align}



The model is fit using basis-specific tuning parameter,
\begin{align}
\widetilde w^{\widehat S, T}_k &= \bigg|\sum_{i=1}^n  R_{ik}^{\widehat S} \frac{\partial \widehat \epsilon_i}{\partial T_i}\bigg|\\
&= \bigg|  R_{k}^{\widehat S, \top}\Delta R^{TX}(c^T-\widehat c^T)\bigg|
\end{align}
since $\partial \epsilon_i/\partial T_i=0$.   Since each element of $\Delta R^{TX}$ is also in $R^X$, by construction, we get the weights to be of order $\sqrt{n \log(p)}$, from which the oracle result follows directly.



\section{Data Appendix}
\label{app:simdata}




\subsection{Larreguy and Marshall Data}

\paragraph{Data Preparation and Transformation}

The original analysis selected a set of state $\times$ covariate interactions to place into a TSLS model.  Rather than select the interactions a priori, we instead create a model of fully saturated interactions and then use the SIS screen to winnow them down.  As the data have a hiearchical structure, we include each variable three times: first, the original variable; second, the mean of the variable by state; and third, the state-centered version of this variable.  Specifically, assume $X_{var}$ are the matrix of variables and $X_{fe}$ the matrix of state fixed effects.  Denote as $H_{fe}$ the hat matrix from the fixed effects and $I$ the commensurate $n \times n$ identity matrix.

The data we place into our software is then
\begin{align}
X^{big} = [X_{fe}:X_{var}: H_{fe} X_{var}:(I-H_{fe}) X_{var}].
\end{align}
All two way interactions and spline transformations of this data are then calculated as normal.

\paragraph{Types of coefficients}

Figure~\ref{fig:coeftype} plots estimates for several different types of variables that get included in our model using data from \citet{larreguy2017}. The top left plots all coefficients that came out of the Sure Independence Screen. The top right plots coefficients on the linear terms, the bottom left on terms that had an interaction, but no non-linear basis function. Finally, the bottom right plots coefficients on variables constructed via an interaction term with a non-linear basis function. We provide Figure~\ref{fig:coeftype} for descriptive purposes. Analysts could use this information, for example, to explore sources of heterogeneity in the LICE.

\begin{figure}[t!]
  \begin{center}\spacingset{1}
    \includegraphics[scale=.5]{CoefficientTypes.pdf}
    \vspace{-.25in}
  \end{center}\spacingset{1}
  \caption{{\bf Coefficient Estimates by Effect Type}      } \label{fig:coeftype}
\end{figure}

\end{document}

mde_arxiv.bbl
\begin{thebibliography}{xx}

\harvarditem{Abadie}{2003}{Abadie2003}
Abadie, Alberto. 2003.
\newblock ``Semiparametric instrumental variable estimation of treatment
  response models.'' {\em Journal of Econometrics} 113(231--263).

\harvarditem[Acharya, Blackwell \harvardand\ Sen]{Acharya, Blackwell
  \harvardand\ Sen}{2016}{acharya2015explaining}
Acharya, Avidit, Matthew Blackwell \harvardand\ Maya Sen. 2016.
\newblock ``Explaining Causal Findings Without Bias: Detecting and Assessing
  Direct Effects.'' {\em American Political Science Review} 110(3).

\harvarditem[Alhamzawi, Yu \harvardand\ Benoit]{Alhamzawi, Yu \harvardand\
  Benoit}{2012}{Alhamzawi2012}
Alhamzawi, Rahim, Keming Yu \harvardand\ Dries~F Benoit. 2012.
\newblock ``Bayesian adaptive Lasso quantile regression.'' {\em Statistical
  Modelling} 12(3):279--297.

\harvarditem[Angrist, Imbens \harvardand\ Rubin]{Angrist, Imbens \harvardand\
  Rubin}{1996}{angrist1996identification}
Angrist, Joshua~D, Guido~W Imbens \harvardand\ Donald~B Rubin. 1996.
\newblock ``Identification of causal effects using instrumental variables.''
  {\em Journal of the American statistical Association} 91(434):444--455.

\harvarditem{Angrist \harvardand\ Pischke}{2009}{angr:pisc:09}
Angrist, Joshua~D. \harvardand\ J\"{o}rn-Steffen Pischke. 2009.
\newblock {\em Mostly Harmless Econometrics: An Empiricist's Companion}.
\newblock Princeton:  Princeton University Press.

\harvarditem{Aronow \harvardand\ Carnegie}{2013}{aronow2013beyond}
Aronow, Peter~M \harvardand\ Allison Carnegie. 2013.
\newblock ``Beyond LATE: Estimation of the average treatment effect with an
  instrumental variable.'' {\em Political Analysis} 21(4):492--506.

\harvarditem{Athey \harvardand\ Imbens.}{2016}{Athey16}
Athey, Susan \harvardand\ Guido Imbens. 2016.
\newblock ``Recursive partitioning for heterogeneous causal effects.'' {\em
  Proceedings of the National Academy of Sciences} 113(27):7353--7360.

\harvarditem[Athey, Tibshirani \harvardand\ Wager]{Athey, Tibshirani
  \harvardand\ Wager}{2016}{athey2016solving}
Athey, Susan, Julie Tibshirani \harvardand\ Stefan Wager. 2016.
\newblock ``Solving Heterogeneous Estimating Equations with Gradient Forests.''
  {\em arXiv preprint arXiv:1610.01271} .

\harvarditem{Austin}{2012}{Austin2012}
Austin, Peter~C. 2012.
\newblock ``Using Ensemble-Based Methods for Directly Estimating Causal
  Effects: An Investigation of Tree-Based G-Computation.'' {\em Multivariate
  Behavioral Research} 47:115--135.

\harvarditem[Belloni et~al.]{Belloni, Chen, Chernozhukov \harvardand\
  Hansen}{2012}{belloni2012sparse}
Belloni, Alexandre, Daniel Chen, Victor Chernozhukov \harvardand\ Christian
  Hansen. 2012.
\newblock ``Sparse models and methods for optimal instruments with an
  application to eminent domain.'' {\em Econometrica} 80(6):2369--2429.

\harvarditem[Bhattacharya et~al.]{Bhattacharya, Pati, Pillai \harvardand\
  Dunson}{2015}{bhattacharya2015dirichlet}
Bhattacharya, Anirban, Debdeep Pati, Natesh~S Pillai \harvardand\ David~B
  Dunson. 2015.
\newblock ``Dirichlet--Laplace priors for optimal shrinkage.'' {\em Journal of
  the American Statistical Association} 110(512):1479--1490.

\harvarditem{Bleakley}{2010}{bleakley2010malaria}
Bleakley, Hoyt. 2010.
\newblock ``Malaria eradication in the Americas: A retrospective analysis of
  childhood exposure.'' {\em American Economic Journal: Applied Economics}
  2(2):1--45.

\harvarditem[Bound, Jaeger \harvardand\ Baker]{Bound, Jaeger \harvardand\
  Baker}{1995}{bound1995problems}
Bound, John, David~A Jaeger \harvardand\ Regina~M Baker. 1995.
\newblock ``Problems with instrumental variables estimation when the
  correlation between the instruments and the endogenous explanatory variable
  is weak.'' {\em Journal of the American statistical association}
  90(430):443--450.

\harvarditem{Buhlmann \harvardand\ Yu}{2002}{Buhlmann02}
Buhlmann, Peter \harvardand\ Bin Yu. 2002.
\newblock ``Analyzing Bagging.'' {\em Annals of Statistics} 30(4):926--961.

\harvarditem{Buhlmann \harvardand\ van~de Geer}{2013}{Buhlmann2013}
Buhlmann, Peter \harvardand\ Sara van~de Geer. 2013.
\newblock {\em Statistics for High-Dimensional Data}.
\newblock Berlin:  Springer.

\harvarditem[Carvalho, Polson \harvardand\ Scott]{Carvalho, Polson \harvardand\
  Scott}{2010}{cps2010horseshoe}
Carvalho, C, N Polson \harvardand\ J Scott. 2010.
\newblock ``The Horseshoe Estimator for Sparse Signals.'' {\em Biometrika}
  97:465--480.

\harvarditem[Chernozhukov, Hansen \harvardand\ Spindler]{Chernozhukov, Hansen
  \harvardand\ Spindler}{2015}{chernozhukovinstrumental}
Chernozhukov, Victor, Christian Hansen \harvardand\ Martin Spindler. 2015.
\newblock ``Instrumental Variables Estimation with Very Many Instruments and
  Controls.''.

\harvarditem[Chipman, George \harvardand\ McCulloch]{Chipman, George
  \harvardand\ McCulloch}{2010}{chipman2010bart}
Chipman, Hugh~A, Edward~I George \harvardand\ Robert~E McCulloch. 2010.
\newblock ``BART: Bayesian additive regression trees.'' {\em The Annals of
  Applied Statistics} pp.~266--298.

\harvarditem[Currie, Durban \harvardand\ Eilers]{Currie, Durban \harvardand\
  Eilers}{2006}{Currie06}
Currie, I.~D., M. Durban \harvardand\ P.~H.~C. Eilers. 2006.
\newblock ``Generalized linear array models with applications to
  multidimensional smoothing.'' {\em Journal of the Royal Statistcal Society,
  Series B} 68(2).

\harvarditem{Davidson \harvardand\ Flachaire}{2008}{Davidson08}
Davidson, Russel \harvardand\ Emmanuel Flachaire. 2008.
\newblock ``The wild bootstrap, tamed at last.'' {\em Journal of Econometrics}
  146(1):162--169.

\harvarditem{de~Boor}{1978}{deBoor78}
de~Boor, C. 1978.
\newblock {\em A Practical Guide to Splines}.
\newblock New York: Springer.

\harvarditem{Dehejia \harvardand\ Wahba}{1999}{dehejia1999causal}
Dehejia, Rajeev~H \harvardand\ Sadek Wahba. 1999.
\newblock ``Causal effects in nonexperimental studies: Reevaluating the
  evaluation of training programs.'' {\em Journal of the American statistical
  Association} 94(448):1053--1062.

\harvarditem{Diamond \harvardand\ Sekhon}{2012}{Diamond2012}
Diamond, Alexis \harvardand\ Jasjeet Sekhon. 2012.
\newblock ``Genetic Matching for Estimating Causal Effects.'' {\em Review of
  Econonics and Statistics} .

\harvarditem{Duflo}{2001}{duflo2001}
Duflo, Esther. 2001.
\newblock ``Schooling and Labor Market Consequences of School Construction in
  Indonesia: Evidence from an Unusual Policy Experiment.'' {\em American
  Economic Review} 91(4):795--813.
\newline\harvardurl{http://www.aeaweb.org/articles?id=10.1257/aer.91.4.795}

\harvarditem{Eilers \harvardand\ Marx}{1996}{Eilers96}
Eilers, Paul H.~C. \harvardand\ Brian~D. Marx. 1996.
\newblock ``Flexible smoothing with B-splines and penalties.'' {\em Statistical
  Science} 11(2):89--121.

\harvarditem{Fan \harvardand\ Lv}{2008}{Fan:2008}
Fan, Jianqing \harvardand\ Jinchi Lv. 2008.
\newblock ``Sure independence screening for ultrahigh dimensional feature
  space.'' {\em Journal of the Royal Statistical Society: Series B}
  70:849--911.

\harvarditem[Fan, Feng \harvardand\ Song]{Fan, Feng \harvardand\
  Song}{2012}{fan2012nonparametric}
Fan, Jianqing, Yang Feng \harvardand\ Rui Song. 2012.
\newblock ``Nonparametric independence screening in sparse
  ultra-high-dimensional additive models.'' {\em Journal of the American
  Statistical Association} .

\harvarditem[Fan, Ma \harvardand\ Dai]{Fan, Ma \harvardand\
  Dai}{2014}{fan2014nonparametric}
Fan, Jianqing, Yunbei Ma \harvardand\ Wei Dai. 2014.
\newblock ``Nonparametric independence screening in sparse
  ultra-high-dimensional varying coefficient models.'' {\em Journal of the
  American Statistical Association} 109(507):1270--1284.

\harvarditem{Friedman}{1991}{friedman1991multivariate}
Friedman, Jerome~H. 1991.
\newblock ``Multivariate adaptive regression splines.'' {\em The Annals of
  Statistics} pp.~1--67.

\harvarditem{Griffin \harvardand\ Brown}{2010}{Griffin2010}
Griffin, J.~E. \harvardand\ P.~J. Brown. 2010.
\newblock ``Inference with normal-gamma prior distributions in regression
  problems.'' {\em Bayesian Analysis} 5(1):171--188.

\harvarditem{Griffin \harvardand\ Brown}{2012}{Griffin2012}
Griffin, J.~E. \harvardand\ P.~J. Brown. 2012.
\newblock ``Structuring shrinkage: some correlated priors for regression.''
  {\em Biometrika} 99(2):481--487.

\harvarditem{Gruber \harvardand\ van~der Laan}{2011}{gruber2011tmle}
Gruber, Susan \harvardand\ Mark~J van~der Laan. 2011.
\newblock ``tmle: An R package for targeted maximum likelihood estimation.''.

\harvarditem{Gu}{2002}{Gu2002}
Gu, Chong. 2002.
\newblock {\em Smoothing spline ANOVA models}.
\newblock Springer series in statistics Springer.

\harvarditem[Gyorfi et~al.]{Gyorfi, Koholor, Krzyzak \harvardand\
  Walk}{2002}{Gyorfi2002}
Gyorfi, Laszlo, Michael Koholor, Adam Krzyzak \harvardand\ Harro Walk. 2002.
\newblock {\em A Distribution-Free Theory of Nonparametric Regression}.
\newblock New York: Springer.

\harvarditem{Haavelmo}{1943}{Haavelmo43}
Haavelmo, Trygve. 1943.
\newblock ``The statistical implications of a system of simultaneous
  equations.'' {\em Econometrica} 11:1--12.

\harvarditem{Hainmueller \harvardand\ Hazlett}{2013}{Hainmueller13}
Hainmueller, Jens \harvardand\ Chad Hazlett. 2013.
\newblock ``Kernel Regularized Least Squares: Reducing Misspecification Bias
  with a Flexible and Interpretable Machine Learning Approach.'' {\em Political
  Analysis} 22(2):143--168.
\newline\harvardurl{+ http://dx.doi.org/10.1093/pan/mpt019}

\harvarditem{Hansen}{2017}{hanseneconometrics}
Hansen, Bruce~E. 2017.
\newblock ``Econometrics.'' Unpublished manuscript.

\harvarditem{Hardle \harvardand\ Stoker}{1989}{Hardle89}
Hardle, Wolfgang \harvardand\ Thomas~M. Stoker. 1989.
\newblock ``Investigating Smooth Multiple Regression by the Method of Average
  Derivatives.'' {\em Journal of American Statistical Association} 84:986--95.

\harvarditem[Hartford et~al.]{Hartford, Lewis, Leyton-Brown \harvardand\
  Taddy}{2016}{taddydeep}
Hartford, Jason, Greg Lewis, Kevin Leyton-Brown \harvardand\ Matt Taddy. 2016.
\newblock ``Counterfactual Prediction with Deep Instrumental Variables
  Networks.'' {\em https://arxiv.org/abs/1612.09596} .

\harvarditem{Heckman \harvardand\ Vytlacil}{1999}{Heckman99}
Heckman, James \harvardand\ Edward Vytlacil. 1999.
\newblock ``Local instrumental variables and latent variable models for
  identifying and bounding treatment effects.'' {\em Proceedings of the
  National Academy of Sciences} 96:4730--4734.

\harvarditem{Heckman \harvardand\ Vytlacil}{2005}{Heckman05}
Heckman, James \harvardand\ Edward Vytlacil. 2005.
\newblock ``Structural Equations, Treatment Effects, and Econometric Policy
  Evaluation.'' {\em Econometrica} 73(3):669--738.

\harvarditem{Heckman \harvardand\ Vytlacil}{2007{\em a}}{Heckman07a}
Heckman, James \harvardand\ Edward Vytlacil. 2007{\em a}.
\newblock ``Econometric Evaluation of Social Programs, Part 1: Causal
  Models,Structural Models, and Econometric Policy Evaluation.'' {\em Handbook
  of Econometrics} 6b:4779--4874.

\harvarditem{Heckman \harvardand\ Vytlacil}{2007{\em b}}{Heckman07b}
Heckman, James \harvardand\ Edward Vytlacil. 2007{\em b}.
\newblock ``Econometric Evaluation of Social Programs, Part II: Using the
  Marginal Treatment Effect to Organize Alternative Econometric Estimators to
  Evaluate Social Programs, and to Forecast their Effects in New
  Environments.'' {\em Handbook of Econometrics} 6b:4875--5143.

\harvarditem[Hill, Weiss \harvardand\ Zhai]{Hill, Weiss \harvardand\
  Zhai}{2011}{Hill2011}
Hill, Jennifer, Christopher Weiss \harvardand\ Fuhua Zhai. 2011.
\newblock ``Challenges With Propensity Score Strategies in a High-Dimensional
  Setting and a Potential Alternative.'' {\em Multivariate Behavioral Research}
  46(3):477--513.

\harvarditem{Hirano \harvardand\ Imbens}{2005}{Hirano2005}
Hirano, Keisuke \harvardand\ Guido Imbens. 2005.
\newblock {\em Applied Bayesian Modeling and Causal Inference from
  Incomplete-Data Perspectives}.
\newblock John Wiley and Sons, Ltd, Chichester, UK chapter The Propensity Score
  with Continuous Treatments.

\harvarditem[Ho et~al.]{Ho, Imai, King \harvardand\
  Stuart}{2007}{ho:imai:king:stua:07}
Ho, Daniel~E., Kosuke Imai, Gary King \harvardand\ Elizabeth~A. Stuart. 2007.
\newblock ``Matching as Nonparametric Preprocessing for Reducing Model
  Dependence in Parametric Causal Inference.'' {\em Political Analysis}
  15(3):199--236.

\harvarditem{Holland}{1986}{holl:86}
Holland, Paul~W. 1986.
\newblock ``Statistics and Causal Inference (with Discussion).'' {\em Journal
  of the American Statistical Association} 81:945--960.

\harvarditem{Horowitz}{2014}{Horowitz14}
Horowitz, Joel. 2014.
\newblock ``Ill-posed inverse problems in economics.'' {\em Annual Review of
  Economics} 6.

\harvarditem{Imai \harvardand\ Van~Dyk}{2012}{imai2012causal}
Imai, Kosuke \harvardand\ David~A Van~Dyk. 2012.
\newblock ``Causal inference with general treatment regimes.'' {\em Journal of
  the American Statistical Association} .

\harvarditem[Imai, Keele \harvardand\ Tingley]{Imai, Keele \harvardand\
  Tingley}{2010}{imai2010general}
Imai, Kosuke, Luke Keele \harvardand\ Dustin Tingley. 2010.
\newblock ``A general approach to causal mediation analysis.'' {\em
  Psychological methods} 15(4):309.

\harvarditem[Imai, Keele \harvardand\ Yamamoto]{Imai, Keele \harvardand\
  Yamamoto}{2010}{imai2010identification}
Imai, Kosuke, Luke Keele \harvardand\ Teppei Yamamoto. 2010.
\newblock ``Identification, inference and sensitivity analysis for causal
  mediation effects.'' {\em Statistical Science} pp.~51--71.

\harvarditem{Imai \harvardand\ Ratkovic}{2014}{CBPS14}
Imai, Kosuke \harvardand\ Marc Ratkovic. 2014.
\newblock ``Covariate Balancing Propensity Score.'' {\em Journal of the Royal
  Statistical Society Series B} 76(1):243--263.

\harvarditem{Imbens \harvardand\ Rubin}{2015}{Imbens2015}
Imbens, Guido~W. \harvardand\ Donald~B. Rubin. 2015.
\newblock {\em Causal inference for statistics, social, and biometical
  sciences}.
\newblock {Cambridge University Press}.

\harvarditem{Imbens \harvardand\ Rosenbaum}{2005}{imbens2005robust}
Imbens, G.W. \harvardand\ P.R. Rosenbaum. 2005.
\newblock ``Robust, accurate confidence intervals with a weak instrument:
  quarter of birth and education.'' {\em Journal of the Royal Statistical
  Society: Series A (Statistics in Society)} 168(1):109--126.

\harvarditem{Kang \harvardand\ Guo}{2009}{Kang2009}
Kang, Jian \harvardand\ Jian Guo. 2009.
\newblock ``Self-adaptive Lasso and its Bayesian Estimation.'' Working Paper.

\harvarditem{Kang \harvardand\ Schafer}{2007}{kang2007demystifying}
Kang, Joseph~DY \harvardand\ Joseph~L Schafer. 2007.
\newblock ``Demystifying double robustness: A comparison of alternative
  strategies for estimating a population mean from incomplete data.'' {\em
  Statistical science} pp.~523--539.

\harvarditem[Kennedy, Lorch \harvardand\ Small]{Kennedy, Lorch \harvardand\
  Small}{2017}{kennedyiv}
Kennedy, Edward~H., Scott~A. Lorch \harvardand\ Dylan~S. Small. 2017.
\newblock ``Robust causal inference with continuous instruments using the local
  instrumental variable curve.'' {\em https://arxiv.org/abs/1607.02566} .

\harvarditem{King \harvardand\ Zeng}{2006}{king2006dangers}
King, Gary \harvardand\ Langche Zeng. 2006.
\newblock ``The dangers of extreme counterfactuals.'' {\em Political Analysis}
  14(2):131--159.

\harvarditem{Kleibergen}{2002}{kleibergen2002pivotal}
Kleibergen, Frank. 2002.
\newblock ``Pivotal statistics for testing structural parameters in
  instrumental variables regression.'' {\em Econometrica} 70(5):1781--1803.

\harvarditem{LaLonde}{1986}{lalonde1986evaluating}
LaLonde, Robert~J. 1986.
\newblock ``Evaluating the econometric evaluations of training programs with
  experimental data.'' {\em The American economic review} pp.~604--620.

\harvarditem{Lam}{2013}{Lam13}
Lam, Patrick. 2013.
\newblock Estimating Individual Causal Effects PhD thesis Harvard University.

\harvarditem{Larreguy \harvardand\ Marshall}{2017}{larreguy2017}
Larreguy, Horacio \harvardand\ John Marshall. 2017.
\newblock ``The Effect of Education on Civic and Political Engagement in
  Non-Consolidated Democracies: Evidence from Nigeria.'' {\em Review of
  Economics and Statistics} .

\harvarditem[Leng, Tran \harvardand\ Nott]{Leng, Tran \harvardand\
  Nott}{2014}{Leng2014}
Leng, Chenlei, Minh-Ngoc Tran \harvardand\ David Nott. 2014.
\newblock ``Bayesian Adaptive LASSO.'' {\em Annals of the Institute of
  Statistical Mathematics} 66(2):221--244.

\harvarditem[MacKinnon et~al.]{MacKinnon, Lockwood, Hoffman, West \harvardand\
  Sheets}{2002}{mackinnon2002comparison}
MacKinnon, David~P, Chondra~M Lockwood, Jeanne~M Hoffman, Stephen~G West
  \harvardand\ Virgil Sheets. 2002.
\newblock ``A comparison of methods to test mediation and other intervening
  variable effects.'' {\em Psychological methods} 7(1):83.

\harvarditem{Meinshausen}{2007}{Meinshausen07}
Meinshausen, Nicolai. 2007.
\newblock ``Relaxed LASSO.'' {\em Computational Statistics and Data Analysis}
  52(1):374--393.

\harvarditem{Morgan \harvardand\ Rubin}{2012}{Morgan12}
Morgan, Kari~Lock \harvardand\ Donald~B. Rubin. 2012.
\newblock ``Rerandomization to improve covariate balance in experiments.'' {\em
  Annals of Statistics} 40(2):1263--1282.

\harvarditem{Newey}{1994}{Newey94}
Newey, Whitney. 1994.
\newblock ``Kernel Estimation of Partial Means and a General Variance
  Estimator.'' {\em Econometric Theory} 10(2):233--253.

\harvarditem{Newey \harvardand\ Powell}{2003}{Newey03}
Newey, Whitney \harvardand\ James Powell. 2003.
\newblock ``Instrumental Variable Estimation of Nonparametric Models.'' {\em
  Econometrica} 71(5):1565--78.

\harvarditem{Newey}{2013}{newey2013nonparametric}
Newey, Whitney~K. 2013.
\newblock ``Nonparametric instrumental variables estimation.'' {\em The
  American Economic Review} 103(3):550--556.

\harvarditem{O'Sullivan}{1986}{OSullivan86}
O'Sullivan, Finbarr. 1986.
\newblock ``A Statistical Perspective on Ill-Posed Inverse Problems.'' {\em
  Statistical Science} 1(4).

\harvarditem{Park \harvardand\ Casella}{2008}{park2008bayesian}
Park, Trevor \harvardand\ George Casella. 2008.
\newblock ``The bayesian lasso.'' {\em Journal of the American Statistical
  Association} 103(482):681--686.

\harvarditem{Pearl}{1995}{pearl1995causal}
Pearl, Judea. 1995.
\newblock ``Causal diagrams for empirical research.'' {\em Biometrika}
  82(4):669--688.

\harvarditem{Pearl}{2000{\em a}}{pearl2009causality}
Pearl, Judea. 2000{\em a}.
\newblock {\em Causality}.
\newblock Cambridge university press.

\harvarditem{Pearl}{2000{\em b}}{Pearl2000}
Pearl, Judea. 2000{\em b}.
\newblock {\em Causality: Models, Reasoning, and Inference}.
\newblock Cambridge: Cambrdige University Press.

\harvarditem{Pearl}{2014}{Pearl2014}
Pearl, Judea. 2014.
\newblock ``Trygve Haavelmo and the Emergence of Causal Calculus.'' {\em
  Econometric Theory} .

\harvarditem{Polley \harvardand\ van~der Laan}{N.d.}{polleysuperlearner}
Polley, Eric \harvardand\ Mark van~der Laan. N.d.
\newblock ``SuperLearner: super learner prediction, 2012.'' {\em URL
  http://CRAN. R-project. org/package= SuperLearner. R package version}.
  Forthcoming.

\harvarditem[Polson, Scott \harvardand\ Windle]{Polson, Scott \harvardand\
  Windle}{2014}{polson2014bayesian}
Polson, Nicholas~G, James~G Scott \harvardand\ Jesse Windle. 2014.
\newblock ``The bayesian bridge.'' {\em Journal of the Royal Statistical
  Society: Series B (Statistical Methodology)} 76(4):713--733.

\harvarditem{Polson \harvardand\ Scott}{2012}{Polson12}
Polson, Nicholas \harvardand\ James Scott. 2012.
\newblock ``Local shrinkage rules, Levy processes and regularized regression.''
  {\em Journal of the Royal Statistical Society, {S}eries {B}} 74(2):287--311.

\harvarditem{Ratkovic \harvardand\ Tingley}{2017}{RatkovicTingley17}
Ratkovic, Marc \harvardand\ Dustin Tingley. 2017.
\newblock ``Sparse Estimation and Uncertainty with Application to Subgroup
  Analysis.'' {\em Political Analysis} 1(25):1--40.

\harvarditem[Ravikumar et~al.]{Ravikumar, Lafferty, Liu \harvardand\
  Wasserman}{2009}{Ravikumar09}
Ravikumar, Pradeep, John Lafferty, Han Liu \harvardand\ Larry Wasserman. 2009.
\newblock ``Sparse Additive Models.'' {\em Journal of the Royal Statistcal
  Society, Series B} 71(5):1009--1030.

\harvarditem{Ridgeway}{1999}{Ridgeway99}
Ridgeway, Greg. 1999.
\newblock ``The state of boosting.'' {\em Computing Science and Statistics}
  31:172--181.

\harvarditem{Robins}{1994}{Robins94}
Robins, James~M. 1994.
\newblock ``Correcting for non-compliance in randomized trials using structural
  nested mean models.'' {\em Communs Statist. Theory Meth.} 23:2379--2412.

\harvarditem[Robins, Rotnitzky \harvardand\ Zhao]{Robins, Rotnitzky
  \harvardand\ Zhao}{1994}{robi:rotn:zhao:94}
Robins, James~M., Andrea Rotnitzky \harvardand\ Lue~Ping Zhao. 1994.
\newblock ``Estimation of Regression Coefficients When Some Regressors are Not
  Always Observed.'' {\em Journal of the American Statistical Association}
  89(427):846--866.

\harvarditem{Robins \harvardand\ Ritov}{1997}{robins1997toward}
Robins, James~M \harvardand\ Ya'acov Ritov. 1997.
\newblock ``Toward a Curse of Dimensionality Appropriate (CODA) Asymptotic
  Theory for Semi-Parametric Models.'' {\em Statistics in Medicine}
  16(3):285--319.

\harvarditem{Rockova \harvardand\ George}{Forthcoming}{Rockova17}
Rockova, Veronica \harvardand\ Edward George. Forthcoming.
\newblock ``The Spike-and-Slab LASSO.'' {\em Journal of American Statistical
  Association} .

\harvarditem{Rosenbaum \harvardand\ Rubin}{1983}{rose:rubi:83}
Rosenbaum, Paul~R. \harvardand\ Donald~B. Rubin. 1983.
\newblock ``The Central Role of the Propensity Score in Observational studies
  for Causal Effects.'' {\em Biometrika} 70(1):41--55.

\harvarditem{Rosenbaum \harvardand\ Rubin}{1984}{rose:rubi:84}
Rosenbaum, Paul~R. \harvardand\ Donald~B. Rubin. 1984.
\newblock ``Reducing Bias in Observational Studies Using Subclassification on
  the Propensity Score.'' {\em Journal of the American Statistical Association}
  79(387):516--524.

\harvarditem{Rubin}{1974}{rubin1974estimating}
Rubin, Donald~B. 1974.
\newblock ``Estimating causal effects of treatments in randomized and
  nonrandomized studies.'' {\em Journal of educational Psychology} 66(5):688.

\harvarditem{Rubin}{2005}{Rubin2005}
Rubin, Donald~B. 2005.
\newblock ``Causal Inference Using Potential Outcomes: Design, Modeling,
  Decisions.'' {\em Journal of the Americal Statistical Association}
  100(469):322--331.

\harvarditem{Smith \harvardand\ Todd}{2005}{smith2005does}
Smith, Jeffrey~A \harvardand\ Petra~E Todd. 2005.
\newblock ``Does matching overcome LaLonde's critique of nonexperimental
  estimators?'' {\em Journal of Econometrics} 125(1):305--353.

\harvarditem{Staiger \harvardand\ Stock}{1994}{staiger1994instrumental}
Staiger, Douglas~O \harvardand\ James~H Stock. 1994.
\newblock ``Instrumental variables regression with weak instruments.''.

\harvarditem{Stein}{1981}{Stein81}
Stein, Charles. 1981.
\newblock ``Estimation of the mean of a multivariate normal distribution.''
  {\em Annals of Statistics} 9:1135--51.

\harvarditem[Stone et~al.]{Stone, Hansen, Kooperberg \harvardand\
  Truong}{1997}{Stone1997}
Stone, Charles~J., Mark~H. Hansen, Charles Kooperberg \harvardand\ Young~K.
  Truong. 1997.
\newblock ``Polynomial Splines and Their Tensor Products in Extended Linear
  Modeling.'' {\em The Annals of Statistics} 25(4):1371--1470.

\harvarditem{van~der Laan \harvardand\ Rose}{2011}{tmlebook}
van~der Laan, Mark~J. \harvardand\ Sherri Rose. 2011.
\newblock {\em Targeted Learning Causal Inference for Observational and
  Experimental Data}.
\newblock Springer.

\harvarditem{VanderWeele}{2015}{vanderweele2015explanation}
VanderWeele, Tyler. 2015.
\newblock {\em Explanation in causal inference: methods for mediation and
  interaction}.
\newblock Oxford University Press.

\harvarditem{Vansteelandt}{2009}{vansteelandt2009estimating}
Vansteelandt, Stijn. 2009.
\newblock ``Estimating direct effects in cohort and case--control studies.''
  {\em Epidemiology} 20(6):851--860.

\harvarditem{Wood}{2006}{Wood06}
Wood, Simon. 2006.
\newblock ``Low-Rank Scale-Invariant Tensor Product Smooths for Generalized
  Additive Mixed Models.'' {\em Biometrics} 62(4):1025--1036.

\harvarditem{Wood}{2016}{Wood16}
Wood, Simon. 2016.
\newblock ``P-splines with derivative based penalties and tensor product
  smoothing of unevenly distributed data.''
  https://arxiv.org/pdf/1605.02446.pdf.

\harvarditem{Zou}{2006}{zouadaptive}
Zou, Hui. 2006.
\newblock ``The Adaptive Lasso and Its Oracle Properties.'' {\em Journal of the
  American Statistical Association} 101(476):1418--1429.

\end{thebibliography}
