[name=Theorem]theorem [name=Proposition]proposition
[name=Assumption,numberwithin=section]assumption

# Causal Estimation of Exposure Shifts with Neural Networks: Evaluating the
Health Benefits of Stricter Air Quality Standards in the US

Mauricio Tec, Kevin Josey, Oladimeji Mudele, Francesca Dominici
Department of Biostatistics
Harvard University
<EMAIL_ADDRESS>Corresponding author.

###### Abstract

In policy research, one of the most critical analytic tasks is to estimate the
causal effect of a policy-relevant shift to the distribution of a continuous
exposure/treatment on an outcome of interest. We call this problem _shift-
response function_ (SRF) estimation. Existing neural network methods involving
robust causal-effect estimators lack theoretical guarantees and practical
implementations for SRF estimation. Motivated by a key policy-relevant
question in public health, we develop a neural network method and its
theoretical underpinnings to estimate SRFs with robustness and efficiency
guarantees. We then apply our method to data consisting of 68 million
individuals and 27 million deaths across the U.S. to estimate the causal
effect from revising the US National Ambient Air Quality Standards (NAAQS) for
$\text{PM}_{2.5}$ from 12 $\mathrm{\upmu g/m^{3}}$ to 9 $\mathrm{\upmu
g/m^{3}}$. This change has been recently proposed by the US Environmental
Protection Agency (EPA). Our goal is to estimate, for the first time, the
reduction in deaths that would result from this anticipated revision using
causal methods for SRFs. Our proposed method, called Targeted Regularization
for Exposure Shifts with Neural Networks (tresnet), contributes to the neural
network literature for causal inference in two ways: first, it proposes a
targeted regularization loss with theoretical properties that ensure double
robustness and achieves asymptotic efficiency specific for SRF estimation;
second, it enables loss functions from the exponential family of distributions
to accommodate non-continuous outcome distributions (such as hospitalization
or mortality counts). We complement our application with benchmark experiments
that demonstrate tresnet’s broad applicability and competitiveness.

##  1 Introduction

The field of causal inference has seen immense progress in the past couple of
decades with the development of _targeted_ doubly-robust methods yielding
desirable theoretical efficiency guarantees on estimates of various causal
effects (Van der Laan et al., 2011; Kennedy, 2016). These advancements have
been recently incorporated into the neural network (NN) literature for causal
inference via _targeted regularization_ (TR) (Shi et al., 2019; Nie et al.,
2021). TR methods produce favorable properties for causal estimation by
incorporating a regularization term into a supervised neural network model.
However, it remains an open task to develop a NN method that specifically
targets the causal effect of a shift in the distribution for a continuous
exposure/treatment variable (Muñoz & Van Der Laan, 2012). We call this problem
_shift-response function_ (SRF) estimation. Many scientific questions can be
formulated as an SRF estimation task (Muñoz & Van Der Laan, 2012). Some
notable examples include estimating the health effects of shifts in the
distribution of environmental, socioeconomic, and behavioral variables (e.g.,
air pollution, income, exercise habits) (Muñoz & Van Der Laan, 2012; Díaz &
Hejazi, 2020; Smith et al., 2023).

Our objective is to develop a neural network technique that addresses a timely
and highly prominent regulatory question. More specifically, the EPA is
currently considering whether or not to revise the National Ambient Air
Quality Standards (NAAQS), potentially lowering the permitted annual-average
$\text{PM}_{2.5}$ concentration from 12 to 11, 10 or 9 $\mathrm{\upmu
g/m^{3}}$. We anticipate that the revision of the NAAQS will ultimately result
in a shift to the distribution of $\text{PM}_{2.5}$ concentrations. Our goal
is to estimate, for the first time, the reduction in deaths that would result
from this anticipated shift using causal methods for SRFs.

Figure 1: Estimated mortality reduction under a _cutoff exposure shift_
lowering the annual $\text{PM}_{2.5}$ in all regions below a given threshold.
Uncertainty bands represent the interquartile range from an ensemble of
networks. _Data source_ : US Medicare claims from 2000–2016.

Contributions We contribute to the public debate informing the US
Environmental Protection Agency (EPA) on the effects of modifying air quality
standards. A preview of the results (fully developed in Section 6) is
presented in Figure 1. The figure presents the estimated reduction in deaths
(%) resulting from various shifts to the distribution of $\text{PM}_{2.5}$
across every ZIP-code in the contiguous US between 2000 and 2016. These shifts
limit the maximum concentration of $\text{PM}_{2.5}$ to the cutoff value for
every ZIP-code that exceeds the cutoff. We vary the cutoff in this SRF between
6 $\mathrm{\upmu g/m^{3}}$ and 16 $\mathrm{\upmu g/m^{3}}$ (x-axis). The
y-axis represents the % reduction in deaths that corresponds with each cutoff
threshold. Notably, a NAAQS threshold of 9 $\mathrm{\upmu g/m^{3}}$, would
have had the effect of decreasing elder mortality by 4%. These findings
present a data-driven perspective on the potential health benefits of the
EPA’s proposal.

In our implementation of this analysis, we present a novel method, called
Targeted Regularization for Exposure Shifts with Neural Networks (tresnet),
which introduces two necessary and generalizable methodological innovations to
the TR literature. First, we use a TR loss that specifically targets SRFs,
ensuring that our estimates retain the properties we have come to expect from
TR methods such as asymptotic efficiency and double robustness (Kennedy,
2016). Given standard regularity conditions, these guarantees imply that the
SRF is consistently estimated when either the outcome model or the density-
ratio model for the exposure shift is correctly specified, and achieves the
best possible efficiency rate when both models are correctly specified.
Second, tresnet accommodates non-continuous outcomes belonging to the
exponential family of distributions (such as mortality counts) that frequently
arise in real-world scenarios, including our motivating application. In
addition to its suitability for our application, we assess the performance of
tresnet in a simulation study tailored for SRF estimation, demonstrating
improvements over neural network methods not specifically designed for SRFs.

Related work Recent papers have most often estimated causal effects relating
air pollution to elder mortality using exposure-response functions (ERFs); see
for example Wu et al. (2020); Bahadori et al. (2022); Josey et al. (2023).
However, none of the methods implemented in these studies target an SRF
estimand. We elaborate on the distinction between ERFs and SRFs in Section 2,
emphasizing the latter’s importance for our motivating application and
informing policy.

Neural network-specific methods for causal inference are divided between works
that aim to estimate individualized effects (e.g., Bica et al. (2020); Yoon et
al. (2018)) and those targeting marginal effects from a population with
efficiency guarantees derived from a specific causal estimand. Several methods
in the latter category—including this work—are based on deriving estimating
equations using a functional quantity unique to the target estimand known as
the _efficient influence function_ (EIF). EIFs have been widely studied for
deriving asymptotically efficient estimators with doubly robust properties
(Kennedy, 2016; Bang & Robins, 2005; Robins, 2000; Bickel et al., 1993). EIFs
are also referred to as Neyman orthogonal scores in the double machine
learning literature (Kennedy, 2022). _Targeted regularization_ (TR) (Shi et
al., 2019) links EIF estimation methods to neural network architectures and
optimization. Recent uses of TR include the dragonnet (Shi et al., 2019),
which introduced TR for targeting the average treatment effect of a binary
treatment, and the vcnet (Nie et al., 2021), for targeting the exposure-
response function (ERF) of a continuous exposure.

Outside of the the neural network literature, SRFs have been studied in the
_stochastic intervention_ and _modified treatment policy_ frameworks for
causal inference (Muñoz & Van Der Laan, 2012; Díaz et al., 2021). Insights
gained from the advancements in this space have guided the formulation of our
applied research question and our implementation of a neural network
architecture.

##  2 Problem Statement: The Causal Effect of an Exposure Shift

The notation we use follows standard language and concepts of the potential
outcomes framework (Imbens & Rubin, 2015). Let $(A,Y,{\bm{X}})$ denote a unit
from the target population ${\mathbb{P}}$, where $A\in{\mathcal{A}}$ is a
continuous exposure/treatment variable, $Y\in{\mathcal{Y}}$ is the outcome of
interest, and ${\bm{X}}\in{\mathcal{X}}$ are covariates. Suppose we obtain a
sample of $n$ iid observations $\\{{\bm{X}}_{i},A_{i},Y_{i}\\}_{i=1}^{n}$. The
_potential outcomes_ notation $Y^{a}$ represents the outcome corresponding to
the exposure/treatment level $a$. The potential outcomes for unobserved
exposure values are called _counterfactuals_ and _factuals_ otherwise. The
_consistency_ condition requires that the factual outcomes satisfy $Y^{A}=Y$
for the observed treatment $A$. The _no interference_ condition requires that
the exposure of one unit does not cause the potential outcomes of another unit
111 _Consistency_ and _no interference_ encompass the _the stable unit
treatment value assignment_ (SUTVA) assumption, which is a standard structural
assumption used to formalize causal estimation problems in the potential
outcomes framework (Imbens & Rubin, 2015)..

The conditional expectation of the potential outcomes is defined as
$\mu({\bm{x}},a)=\mathbb{E}[Y^{a}|{\bm{X}}={\bm{x}}]$. All expectations are
with respect to ${\mathbb{P}}$ unless stated otherwise. For conciseness, we
also use a generic notation $p$ to denote the density function of the random
variables composing the data distribution. For example, ${\bm{X}}\sim
p({\bm{X}})$, $A|{\bm{X}}\sim p(A|{\bm{X}})$, and so on. The latter quantity
$p(A|{\bm{X}})$ is known as the _generalized propensity score_.

Estimand of interest: effects of exposure shifts An exposure shift represents
a counterfactual scenario in which an intervention modifies a unit’s exposure,
usually in reference to its observed value. It considers the case when the
exposure data consists of pairs $(A,\tilde{A})$ in which $A$ is the exposure
observed in the data and $\tilde{A}$ indicates the exposure that would be
observed after the shift. As simple examples, illustrated in Figure 2, one can
define a cutoff shift $\tilde{A}=\min\\{A,c\\}$ where the exposure is
truncated to a maximum value of $c$; another example is $\tilde{A}=cA$, in
which the exposure is multiplied by a factor $c$. For instance, $c=0.9$ would
represent a $10\%$ reduction to all units with respect to its observed value.
We will denote the shifted generalized propensity score as
$\tilde{p}(\tilde{A}|{\bm{X}})$.

The shift-response function (SRF) estimand, denoted $\psi=\psi({\mathbb{P}})$,
is the expected potential outcome induced by the exposure shift:

$\psi=\mathbb{E}_{{\bm{X}}\sim
p({\bm{X}})}\left[\mathbb{E}_{\tilde{A}\sim{\tilde{p}(\tilde{A}|{\bm{X}})}}\left[{\mu({\bm{X}},\tilde{A})}\mid{\bm{X}}\right]\right]={\mathbb{E}[Y^{\tilde{A}}]}.$
(1)

This estimand cannot be expressed in terms of traditional causal effects such
as the average treatment effect (ATE) or an exposure-response function (ERF)
(Muñoz & Van Der Laan, 2012).

There is no restriction on how $\tilde{A}$ is defined as long as pairs
$(A,\tilde{A})$ are available. The shift can be a stochastic or deterministic
function $\tilde{A}=f(A,{\bm{X}})$ (Díaz et al., 2021). Importantly, we do not
require $\tilde{p}$ to be known, nor do we need a deterministic formula
expressing $\tilde{A}$ in terms of $A$.

(a) Cutoff shift

(b) Percent reduction shift

(c) ERF

Figure 2: Two examples of exposure shifts with their implied counterfactuals
and, for comparison, the implied counterfactuals of an exposure-response
function at a given exposure value.

Causal identification The target estimand $\psi$ can be expressed as a
functional of the observable data distribution under standard assumptions,
which are:

[Unconfoundedness] $A\perp\\!\\!\\!\\!\perp Y^{a}\mid{\bm{X}}$ for all
$a\in{\mathcal{A}}$. [Positivity] Let
$w({\bm{x}},a)=\tilde{p}(a|{\bm{x}})/p(a|{\bm{x}})$. Then, there exist a
constant $M>0$ such that $w({\bm{x}},a)<M$ for all $(a,{\bm{x}})$ such that
$p(a|{\bm{x}})>0$.

The first assumption ensures that
$\mu({\bm{x}},a)=\mathbb{E}[Y|{\bm{X}}={\bm{x}},A=a]$. The right-hand side of
this equality can be directly estimated via regression whereas the left-hand
side cannot. The second assumption implies that the density ratio
$w({\bm{x}},a)$ is well-defined and behaved. Notice that
$\psi=\mathbb{E}[\mu({\bm{X}},A)w({\bm{X}},A)]$ by the importance sampling
formula. Therefore, estimators of $\mu$ and $w$ will suffice to estimate
$\psi$. Intuitively, Figure 2 prohibits extreme cases when the shifted
exposures take value outside the practical domain of the observed exposure, in
which case counterfactual estimation is impossible (Muñoz & Van Der Laan,
2012).

Mutiple shifts We can estimate the effect of multiple exposure shifts
simultaneously. Let $\tilde{p}\in\tilde{{\mathcal{P}}}$ denote the set of
finite exposure shifts of interest. We can index the density ratio and
estimand by $\tilde{p}$, and denote
${\bm{w}}=(w_{\tilde{p}})_{\tilde{p}\in\tilde{{\mathcal{P}}}}$,
${\bm{\psi}}=(\psi_{\tilde{p}})_{\tilde{p}\in\tilde{{\mathcal{P}}}}$ and
$\tilde{{\bm{A}}}=(\tilde{A}^{\tilde{p}})_{\tilde{p}\in\tilde{P}}$.

Comparison with ERFs Exposure-response functions are common estimands in the
causal inference literature for continuous treatments. Mathematically, the ERF
$\xi$ can be written as the mapping $\xi(a)=\mathbb{E}[\mu({\bm{X}},A)|A=a]$.
One can consider ERFs as a limiting case of SRFs when $\tilde{p}$ is a point
mass distribution centered at a fixed treatment value assigned _equally to all
units_. A visual example is shown in Figure 2(c). The fundamental reason why
SRFs are more suitable for our motivating application is that ERFs do not
allow us to consider scenarios like the cutoff shift, illustrated in Figure 1.
In this setting the $\text{PM}_{2.5}$ levels are reduced only for locations
that did not comply with the proposed air quality standard. An ERF describes
the average outcome when _all_ units are given the _same_ $\text{PM}_{2.5}$
value. Thus, theoretical guarantees for ERF estimation do not apply to SRFs.
Moreover, our experiments suggest that the TR estimators designed for ERFs
hampers estimation of SRFs.

##  3 tresnet: Targeted Regularization for Estimating the Causal Effects of
Exposure Shifts with Neural Networks

As we have suggested earlier, an estimator of ${\bm{\psi}}$ can be derived
from estimators of the outcome and density ratio functions. Using TR, we will
obtain an estimator $\hat{{\bm{\psi}}}^{\text{tr}}$ with the guarantees that
$\lVert\hat{{\bm{\psi}}}^{\text{tr}}-{\bm{\psi}}\rVert_{2}$ converges in
probability at an “efficient” rate according to the prevailing semiparametric
efficiency theory surrounding robust causal effect estimation (Kennedy, 2022).

The efficient influence function (EIF) The EIF of ${\bm{\psi}}$, denoted
$\varphi({\bm{O}};{\bm{\psi}},\mu,{\bm{w}})$, is a fundamental function in the
theory of semiparametric models (Kennedy, 2022). More concisely, the EIF is
the gradient of ${\bm{\psi}}$ with respect to small perturbations in the data
distribution. Results from semiparametric theory show that the best possible
variance among the family of regular, asymptotically linear estimators of
${\bm{\psi}}$ is bounded below by
${\mathbb{P}}[\varphi({\bm{\psi}},\mu,{\bm{w}})\varphi({\bm{\psi}},\mu,{\bm{w}})^{\top}]$.
Moreover, the asymptotic variance of any statistically consistent estimator
$(\hat{{\bm{\psi}}},\hat{\mu},\hat{{\bm{w}}})$ satisfying the empirical
estimating equation (EEE)
${\mathbb{P}}_{n}(\varphi(\hat{{\bm{\psi}}},\hat{\mu},\hat{{\bm{w}}}))=0$
eventually achieves this lower bound. For the SRF, the EIF is given by

$\displaystyle\varphi({\bm{O}};{\bm{\psi}},\mu,{\bm{w}})={\bm{w}}({\bm{X}},A)\left(Y-\mu({\bm{X}},A)\right)+\mu({\bm{X}},\tilde{{\bm{A}}})-{\bm{\psi}}.$
(2)

We provide a proof and additional background about the derivation of the EIF
in the appendix. The reader can refer to Tsiatis (2006) and Kennedy (2022) for
a more comprehensive introduction to the EIF and semiparametric efficiency
theory. Observe that if $(\hat{{\bm{\psi}}},\hat{\mu},\hat{{\bm{w}}})$
satisfies the EEE, then $\hat{{\bm{\psi}}}$ can be decomposed in terms of a
debiasing component of the residual error and a plugin estimator for the
marginalized average of the mean response:

$\hat{{\bm{\psi}}}={\color[rgb]{0,0,1}\underbrace{\textstyle{\frac{1}{n}\sum_{i=1}^{n}\hat{{\bm{w}}}({\bm{X}}_{i},A_{i})(Y_{i}-\hat{\mu}({\bm{X}}_{i},A_{i}))}}_{\text{debiasing
term}}}+{\color[rgb]{1,0,0}\underbrace{\textstyle{\frac{1}{n}\sum_{i=1}^{n}\hat{\mu}({\bm{X}}_{i},\tilde{{\bm{A}}}_{i})}}_{\text{\color[rgb]{1,0,0}
plugin estimator}}}$ (3)

TR for SRFs An immediate approach to obtain a doubly-robust estimator
satisfying the EEE would be to use the right-hand side of Equation 3 as
specifying an estimator from a finite sample and nuisance function estimators
$\hat{\mu}$ and $\hat{{\bm{w}}}$. Such an estimator, denoted
$\hat{{\bm{\psi}}}_{\text{aipw}}$, is sometimes called the _augmented inverse-
probability weighting_ (AIPW) estimator for exposure shifts (also referred to
as modified treatment policies and stochastic interventions) (Muñoz & Van Der
Laan, 2012; Díaz et al., 2021). This estimator must be distinguished from the
standard AIPW estimator for traditional average causal effects (ATE and ERF)
(Robins, 2000; Robins et al., 2000).

TR is an alternative approach based on the observation that the debiasing term
in Equation 3 has been empirically observed to affect performance in finite
samples due to its sensitivity to $\hat{{\bm{w}}}$. Instead, TR learns a
perturbed outcomes model using a special regularization loss, ensuring that
the resulting plugin estimator (the second component of Equation 3) satisfies
the EEE without requiring the debiasing term. We introduce the perturbation
model and regularization loss in the next section.

Generalized TR for outcomes in the exponential family We present a general
formulation applicable to the SRF estimand for any outcome supported by a
generalized domain. First, we say that the outcome follows a conditional
distribution from the exponential family if
$p(Y|{\bm{X}},A)\propto\exp(Y\eta({\bm{X}},A)-\Lambda(\eta({\bm{X}},A))$ for
some function $\eta:{\mathcal{X}}\times{\mathcal{A}}\to\mathbb{R}$. The
family’s canonical link function $g$ is defined by the identity
$g(\mathbb{E}[Y\mid{\bm{X}},A])=\eta({\bm{X}},A)$. For all distributions in
the exponential family, $g$ is invertible (McCullagh, 2019). Exponential
families allow us to consider the usual mean-squared error and logistic
regression as special cases. They also enable modeling of death counts as in
our application analyzing the health effects of $\text{PM}_{2.5}$. In this
setting, we set $\Lambda(\eta)=e^{\eta}$ and $g(\mu)=\log(\mu)$; a Poisson
regression environment.

The following theorem forms the basis for the TR estimator.

[] Let ${\bm{\epsilon}}$ denote a perturbation parameter and define

$\displaystyle{\mathcal{L}}^{\text{tr}}(\mu^{\text{NN}},{\bm{w}}^{\text{NN}},{\bm{\epsilon}})({\bm{O}})$
$\displaystyle=\Lambda(g(\mu^{\text{NN}}({\bm{X}},A))+{\bm{\epsilon}})-(g(\mu^{\text{NN}}({\bm{X}},A))+{\bm{\epsilon}})Y.$
(4)
$\displaystyle{\mathcal{R}}^{\text{tr}}(\mu^{\text{NN}},{\bm{w}}^{\text{NN}},{\bm{\epsilon}})$
$\displaystyle=\textstyle{\frac{1}{n}\sum_{i=1}^{n}{\mathcal{L}}^{\text{tr}}(\mu^{\text{NN}},{\bm{w}}^{\text{NN}},{\bm{\epsilon}})({\bm{O}}_{i})}.$

Then
$(\frac{\partial{\mathcal{R}}^{\text{tr}}}{\partial{\bm{\epsilon}}})({\mu}^{\text{NN}},{{\bm{w}}}^{\text{NN}},{{\bm{\epsilon}}})=0$
iff
$\frac{1}{n}\sum_{i=1}^{n}{{\bm{w}}}^{\text{NN}}({\bm{X}}_{i},A_{i})(Y_{i}-g^{-1}(g({\mu}^{\text{NN}}({\bm{X}}_{i},A_{i}))+{{\bm{\epsilon}}})))=0$.

The condition
${\partial{\mathcal{R}}^{\text{tr}}}/{\partial{\bm{\epsilon}}}=0$ in the
theorem holds upon minimization of ${\mathcal{R}}^{\text{tr}}$. Consequently,
the TR estimator $\hat{{\bm{\psi}}}^{\text{tr}}$ is defined as the solution of
an optimization problem

$\displaystyle(\hat{\mu},\hat{{\bm{w}}},\hat{{\bm{\epsilon}}})$
$\displaystyle=\operatorname*{arg\,min}_{\mu^{\text{NN}},{\bm{w}}^{\text{NN}},{\bm{\epsilon}}}{\mathcal{R}}_{\mu}(\mu^{\text{NN}})+\alpha{\mathcal{R}}_{\bm{w}}({\bm{w}}^{\text{NN}})+\beta_{n}{\mathcal{R}}^{\text{tr}}(\mu^{\text{NN}},{\bm{w}}^{\text{NN}},{\bm{\epsilon}})$
(5) $\displaystyle\hat{{\bm{\psi}}}^{\text{tr}}$
$\displaystyle:=\textstyle{\frac{1}{n}\sum_{i=1}^{n}g^{-1}(g(\hat{\mu}({\bm{X}},A))+\hat{{\bm{\epsilon}}}))}$

where ${\mathcal{R}}_{\mu}$ and ${\mathcal{R}}_{\bm{w}}$ are the empirical
risk functions of $\mu$ and ${\bm{w}}$, $\alpha>0$ is a hyperparameter, and
$\beta_{n}$ is a regularization weight satisfying $\beta_{n}\to 0$. The latter
condition is needed to ensure statistical consistency, as first discussed by
Nie et al. (2021) for ERF estimation. Section 4 provides additional details
about the architecture and risk function specification. The full loss in
Equation 5 preserves the fact that
$\frac{\partial{\mathcal{R}}^{\text{tr}}}{\partial{\bm{\epsilon}}}=0$ since
${{\bm{\epsilon}}}$ only appears in the regularization term.

Double robustness and efficiency Before introducing the main result of the TR
estimator for the SRF, we require the following additional notation. For any
$f:{\mathcal{U}}\to\mathbb{R}$, $\lVert
f\rVert_{\infty}=\sup_{u\in{\mathcal{U}}}|f(u)|$ and
$\lVert{\mathcal{F}}\rVert_{\infty}=\sup_{f\in{\mathcal{F}}}\lVert
f\rVert_{\infty}$. We define the sample Rademacher complexity as
$\mathrm{Rad}_{n}({\mathcal{F}})=\sup_{f\in{\mathcal{F}}}|\frac{1}{n}\sum_{i=1}^{n}\sigma_{i}f(U_{i})|$
where $\sigma_{i}$ are iid Rademacher random variables satisfying
$p(\sigma_{i}=1)=p(\sigma_{i}=-1)=1/2$. It is easiest to think of the
Rademacher complexity as a natural measure for the degrees of freedom of an
estimating class of functions. We use $O_{p}$ and $o_{p}$ to denote stochastic
boundedness and convergence in probability, respectively.

[] Let ${\mathcal{M}}$ and ${\mathcal{W}}$ be classes of functions such that
${\hat{\mu},\mu}\in{\mathcal{M}}$ and
${\hat{{\bm{w}}},{\bm{w}}}\in{\mathcal{W}}$. Suppose assumptions 2 and 2 hold,
and that the following regularity conditions hold: (i)
$\lVert{\mathcal{M}}\rVert_{\infty}<\infty$,
$\lVert{\mathcal{W}}\rVert_{\infty}<\infty$, $\lVert
1/{\mathcal{W}}\rVert_{\infty}<\infty$; (ii) either $\hat{\mu}=\mu$,
$\hat{{\bm{w}}}={\bm{w}}$, or $\mathrm{Rad}_{n}({\mathcal{M}})=O(n^{-1/2})$
and $\mathrm{Rad}_{n}({\mathcal{W}})=O(n^{-1/2})$; (iii) the loss function in
Equation 4 is Lipschitz; (iv) $\Lambda$ and $g$ are twice continuously
differentiable. Then, the following statements are true:

  1. 1.

The outcome and density ratio estimators of TR are consistent. That is,
$\hat{\mu}{\overset{p}{\longrightarrow}}\mu$ and
$\hat{{\bm{w}}}{\overset{p}{\longrightarrow}}{\bm{w}}$.

  2. 2.

The estimator $\hat{{\bm{\psi}}}^{\text{tr}}$ satisfies
$\lVert\hat{{\bm{\psi}}}^{\text{tr}}-{\bm{\psi}}\rVert_{\infty}=O_{p}(n^{-1/2}+r_{1}(n)r_{2}(n))$
whenever $\lVert\hat{\mu}-\mu\rVert_{\infty}={O_{p}}(r_{1}(n))$ and
$\lVert\hat{{\bm{w}}}-{\bm{w}}\rVert_{\infty}={O_{p}}(r_{2}(n))$.

Section 3 shows that the TR regularized learner of the SRF achieves “optimal”
root-$n$ convergence when $r_{1}(n)=r_{2}(n)=n^{-1/4}$ or when either $r_{1}$
or $r_{2}$ vanishes. Using standard arguments involving concentration
inequalities, the Lipschitz assumption on the loss function can be relaxed by
assuming that the loss function has a vanishing Rademacher complexity
(Wainwright, 2019).

Figure 3: tresnet architecture using a head for the density ratio model and a
head for the outcome model.

##  4 Architecture for Estimating $\mu$ and $w$ with Neural Networks

We describe a simple yet effective architecture for estimating $\mu$ and
${\bm{w}}$. To keep the notation simple, we will write $f_{\theta}$ to denote
generic output from a neural network indexed by weights $\theta$. The
architecture has three components, illustrated in Figure 3. The first
component maps the confounders ${\bm{X}}$ to a latent representation
${\bm{Z}}=f_{\theta_{Z}}({\bm{X}})\in\mathbb{R}^{d}$. This component will
typically be a multi-layer perception (MLP). The second and third components
are the outcome and density-ratio heads, which are functions of ${\bm{Z}}$ and
the treatment, respectively. We describe all three components in detail below.

Neural network architectures for nuisance function estimation have been widely
investigated in causal inference; see Farrell et al. (2021) for a review. We
use the architectures proposed in the TR literature as a building block,
particularly for continuous treatments (Nie et al., 2021). Nonetheless, It
must be remarked that previous work in TR has not yet investigated
architectures required for SRF estimation. In particular, we need a new
architecture to estimate the density ratio ${\bm{w}}$. Previous works have
only focused on architectures for estimating propensity scores as required by
traditional causal effect estimation (ATEs and ERFs).

Outcome model Recall that we assume the outcome $Y$ follows a conditional
distribution from the exponential family. That is,
$p(Y|{\bm{X}},A)\propto\exp(Y\eta({\bm{X}},A)-\Lambda(\eta({\bm{X}},A))$ with
an invertible link function $g$ satisfying
$\mu({\bm{X}},A)=g^{-1}(\eta({\bm{X}},A))$. We can identify the canonical
parameter $\eta$ with the output of the neural network and learn $\hat{\mu}$
by minimizing the empirical risk

$\displaystyle{\mathcal{R}}_{\mu}(\mu^{\text{NN}})=\frac{1}{n}\sum_{i=1}^{n}\left\\{\Lambda(g^{-1}(\mu^{\text{NN}}({\bm{Z}},A))-Yg^{-1}(\mu^{\text{NN}}({\bm{Z}},A))\right\\}.$
(6)

Next, we need to select a functional form for the neural network. An MLP
parameterization with the concatenated inputs of $({\bm{Z}},A)$–the naïve
choice–would likely result in the effect of $A$ being lost in intermediate
computations. Instead, we adopt the _varying coefficient_ approach by setting
$\mu^{\text{NN}}({\bm{X}},A)=g^{-1}(f_{\theta_{\mu}(A)}({\bm{Z}}))$ (Nie et
al., 2021; Chiang et al., 2001). With this choice, the weights of each layer
are dynamically computed as a function of ${\bm{A}}$ obtained from a linear
combination of basis functions spanning the set of admissible functions on
${\bm{A}}$. The weights of the linear combination are themselves a learnable
linear combination of the hidden outputs from the previous layer. We refer the
reader to Nie et al. (2021) for additional background on varying-coefficient
layers. Our experiments suggest TR is beneficial for different choices of
basis functions.

### Estimation of ${\bm{w}}$ via classification

The density ratio head ${\bm{w}}^{\text{NN}}$ is trained using an auxiliary
classification task. The goal is to estimate the density ratio $w_{j}$ for
each $j=1,\ldots,|\tilde{P}|$. For this purpose, we use an auxiliary
classification task where the positive labels are assigned to the samples from
$\tilde{A}_{j}$ and the negative labels to the samples with $A$ such that

${\mathcal{R}}_{w}({\bm{w}}^{\text{NN}})=\sum_{i=1}^{n}\sum_{\tilde{p}\in\tilde{P}}\frac{1}{2n|\tilde{P}|}\left\\{\textrm{BCE}(\log
w^{\text{NN}}_{\tilde{p}}({\bm{X}},\tilde{A}_{i\tilde{p}}),{\bm{1}})+\textrm{BCE}(\log
w^{\text{NN}}_{\tilde{p}}({\bm{X}}_{i},A_{i}),{\bm{0}})\right\\}$ (7)

where BCE stands for the binary cross-entropy classification loss. Equation 7
is a multi-head neural network version of the loss proposed in Díaz et al.
(2021) to estimate the effects of modified treatment policies. To capture the
role of $A$ more accurately, we propose to parameterize the network using the
varying-coefficient structure discussed in the previous section with $\log
w_{\tilde{p}}^{\text{NN}}({\bm{X}},A)=f_{\theta^{\tilde{p}}_{w}(A)}({\bm{Z}})$.
To our knowledge, we are the first to consider a varying-coefficient
architecture for density ratio estimation.

| Spline-based varying coefficients | Piecewise linear varying coefficients
---|---|---
Dataset | aipw${}_{\text{vc}}$ | outcome${}_{\text{vc}}$ | tresnet${}_{\text{vc}}$ | vcnet | aipw${}_{\text{pl}}$ | drnet | tresnet${}_{\text{pl}}$ | drnet+tr${}_{\text{erf}}$
ihdp | 3.15 (0.37) | 2.19 (0.06) | 0.61 (0.03) | 0.63 (0.03) | 1.18 (0.14) | 2.36 (0.06) | 0.15 (0.02) | 0.19 (0.02)
news | 1.5 (0.19) | 3.65 (0.04) | 0.18 (0.02) | 0.28 (0.03) | 0.99 (0.12) | 0.99 (0.1) | 0.17 (0.01) | 0.26 (0.03)
sim-B | 4.1 (0.57) | 0.5 (0.05) | 0.26 (0.03) | 0.29 (0.04) | 1.46 (0.2) | 1.6 (0.2) | 0.14 (0.02) | 0.16 (0.02)
sim-N | 5.69 (0.64) | 0.52 (0.05) | 0.32 (0.02) | 0.32 (0.03) | 1.81 (0.25) | 0.95 (0.06) | 0.14 (0.01) | 0.15 (0.01)
tcga-1 | 1.13 (0.08) | 0.63 (0.02) | 0.8 (0.01) | 0.87 (0.03) | 0.76 (0.05) | 0.62 (0.02) | 0.61 (0.02) | 0.69 (0.03)
tcga-2 | 0.76 (0.09) | 0.24 (0.02) | 0.18 (0.01) | 0.24 (0.02) | 0.36 (0.05) | 0.17 (0.01) | 0.12 (0.0) | 0.16 (0.01)
tcga-3 | 0.83 (0.11) | 0.38 (0.03) | 0.1 (0.01) | 0.15 (0.02) | 0.59 (0.06) | 0.59 (0.04) | 0.08 (0.01) | 0.15 (0.02)

(a) Performance of tresnet for two baseline architectures, including
comparisons with vcnet and drnet.

experiment | outcome${}_{\text{vc}}$w/poisson loss | outcome${}_{\text{vc}}$w/mse loss | tresnet${}_{\text{vc}}$w/poisson loss | tresnet${}_{\text{vc}}$w/mse loss
---|---|---|---|---
ihdp | 18.82 (3.18) | 3726.92 (357.31) | 2.04 (0.08) | 10986.43 (211.99)
news | 3.41 (0.25) | 372.24 (52.02) | 0.33 (0.05) | 1187.94 (100.8)
sim-B | 1222.61 (1269.53) | 8433.98 (1327.22) | 1113.58 (1270.49) | 16902.41 (1233.54)
sim-N | 50.72 (4.45) | 2491.63 (310.36) | 4.6 (0.22) | 18062.85 (243.2)
tcga-1 | 184.89 (6.67) | 6682.91 (474.32) | 40.56 (19.2) | 16307.15 (232.64)
tcga-2 | 48.3 (1.44) | 5854.08 (483.41) | 398.82 (143.96) | 16112.02 (186.99)
tcga-3 | 18.27 (2.73) | 14381.06 (516.25) | 12.92 (2.3) | 8565.17 (447.13)

(b) Performance of Poisson-based regularization when the true data is Poisson
count data

Table 1: Experiment results. The table shows the $\sqrt{\textsc{mise}}$ across
100 random seeds with 95% confidence intervals computed with the asymptotic
normal formula.

##  5 Simulation study

We conducted simulation experiments to validate the design choices of tresnet.

Synthetic and semi-synthetic benchmark datasets The so-called fundamental
problem of causal inference is that the counterfactuals are never observed in
real data. Thus, we need to rely on widely used semi-synthetic datasets to
evaluate the validity of our proposed estimators. First, we consider two
datasets introduced by Nie et al. (2021), which are continuous-treatment
adaptions to the popular datasets ihdp (Hill, 2011) and news (Newman, 2008).
We also use Nie et al. (2021)’s fully simulated data, sim-N. In addition, we
consider the fully simulated dataset described in Bahadori et al. (2022),
which features a continuous treatment and has been previously used for
calibrating models in air pollution studies. Finally, we consider three
variants of the tcga dataset presented in Bica et al. (2020). The three
variants consist of three different dosage specifications as the treatment
assignment and the corresponding dose-response as the outcome. The datasets
described here have been employed without substantial modifications from the
original source studies to facilitate fair comparison. Note that for each of
these synthetic and semi-synthetic datasets, we have access to the true
counterfactuals, which allows us to compute SRFs exactly. We will consider the
estimation task of 20 equally-spaced percent reduction shifts between 0-50%
from the current observed exposures. More specifically, $\tilde{A}=(1-c)A$ for
values of $c$ in $0-50\%$.

Evaluation metric and task Given a semi-synthetic dataset ${\mathcal{D}}$,
consider an algorithm that produces an estimator
$\hat{{\bm{\psi}}}^{(s)}_{\tilde{p},{\mathcal{D}}}$ of
${{\bm{\psi}}}^{(s)}_{\tilde{p},{\mathcal{D}}}$ given an exposure shift
$\tilde{p}$ and random seed $s$. To evaluate the quality of the estimator, we
use the _mean integrated squared error_
$\textsc{mise}_{\mathcal{D}}=\frac{1}{n_{\text{seeds}}|\tilde{{\mathcal{P}}}|}\sum_{s}\sum_{\tilde{p}}\lvert\hat{{\bm{\psi}}}^{(s)}_{\tilde{p},{\mathcal{D}}}-{\bm{\psi}}^{(s)}_{\tilde{p},{\mathcal{D}}}\rvert^{2}.$
This metric is a natural adaption of an analogous metric commonly used in
dose-response curve estimation (Bica et al., 2020).

Experiment 1: Does targeting SRFs help? We evaluate two variants of tresnet
against alternative SRF estimators, including vcnet (Nie et al., 2021) and
drnet (Schwab et al., 2020), two prominent methods used in causal TR
estimation. The first variant of tresnet uses varying-coefficient layers based
on splines—see the discussion in Section 4 for background. We compare this
variant, named tresnet${}_{\textsc{vc}}$, with the following baselines: (a)
vcnet, which uses a similar architecture and plugin estimator as
tresnet${}_{\textsc{vc}}$, but with a TR designed for ERFs; (b)
$\textsc{aipw}_{\textsc{vc}}$ which is the doubly robust, augmented inverse
probability weighting (AIPW) estimator for SRFs (Muñoz & Van Der Laan, 2012)
wherein we fit separate outcome and density ratio models that are then
substituted into Equation 3; and (c) $\textsc{outcome}_{\textsc{vc}}$, which
is the same as tresnet${}_{\textsc{vc}}$ but without the TR and density ratio
heads. The second variant of tresnet uses varying coefficients based on
piecewise linear functions instead of splines. This variant,
tresnet${}_{\textsc{pl}}$, is compared against analogous baselines: (a)
$\textsc{aipw}_{\textsc{pl}}$ for the AIPW estimator; (b) drnet, which is
similar to $\textsc{outcome}_{\textsc{vc}}$ with piecewise linear basis
functions; and (c) $\textsc{drnet}+\textsc{tr}_{\textsc{erf}}$, which uses the
same regularization loss for ERFs as vcnet.

Figure 4(b) shows the results of this experiment. For both architectures, the
tresnet variants achieve the best performance. tresnet${}_{\textsc{vc}}$ is
somewhat better than vcnet, which have comparable architectures although
tresnet uses a different TR implementation specific for SRF estimation rather
than ERF estimation. Likewise, tresnet${}_{\textsc{pl}}$ outperforms
$\textsc{drnet}+\textsc{tr}_{\textsc{erf}}$. These moderate but consistent
performance gains suggest the importance of SRF-specific forms of TR. We also
see strong advantages against outcome-based predictions and AIPW estimators,
suggesting that the TR loss and the shared learning architecture is a boon to
performance. These results are compatible with observations from previous work
in the TR literature (Nie et al., 2021; Shi et al., 2019).

Experiment 2: Does TR for count data improve estimation when count-valued
outcomes are observed? We evaluate whether tresnet with the Poisson-specific
TR, explained in Section 3, performs better than the mean-squared error (MSE)
loss variant when the true data follows a Poisson distribution. This
evaluation is important since our application consists of count data and the
Poisson model is widely used to investigate the effects of $\text{PM}_{2.5}$
on health (Wu et al., 2020; Josey et al., 2022). We construct similar semi-
synthetic datasets as in Experiment 1, but in this experiment the outcome-
generating mechanisms samples from a Poisson distribution rather than a
Gaussian distribution. The results of this experiment are clear–using the
correct exponential family for the outcome model is crucial, regardless of
whether TR is implemented. For these experiments, we used the spline-based
variant of tresnet.

##  6 Application: The Effects of Stricter Air Quality Standards

We implemented tresnet for count data to estimate the health benefits caused
by shifts to the distribution of $\text{PM}_{2.5}$ that would result from
lowering the NAAQS—the regulatory threshold for the annual-average
concentration of $\text{PM}_{2.5}$ enforced by the EPA.

Figure 4: Fraction (%) of observed units remaining above $\text{PM}_{2.5}$
limit as a function of reduction (%) considering different NAAQS (current
NAAQS is set at 12 $\mathrm{\upmu g/m^{3}}$).

Data The dataset is comprised of Medicare data222Access to Medicare data is
restricted without authorization by the Centers for Medicare & Medicaid
Services since it contains sensitive health information. The authors have
completed the required IRB and ethical training to handle these datasets. from
2000–2016, involving 68 million individuals. The data includes measurements on
participant race/ethnicity, sex, age, Medicaid eligibility, and date of death,
which are subsequently aggregated to the annual ZIP-code level. The
$\text{PM}_{2.5}$ exposure measurements are extracted from an ensemble
prediction model (Di et al., 2019). The confounders include measurements on
meteorological information, demographics, and the socioeconomic status of each
ZIP-code. Calendar year and census region indicators are also included to
account for temporal and spatial trends. To compile our dataset, we replicated
the steps and variables outlined by Wu et al. (2020).

Exposure shifts We consider two types of $\text{PM}_{2.5}$ shifts, _cutoff
shifts_ and _percent reduction shift_ : each providing different perspectives
and insights. The counterfactuals implied from these scenarios are illustrated
in figures 2(a) and 2(b), respectively. First, a _cutoff shift_ ,
parameterized by a threshold $d$, encapsulates scenarios in which every ZIP-
code year that exceeded some threshold are truncated to that maximum
threshold. Mathematically, the shift is defined by the transformation
$\tilde{A}=\min(A,c)$. To be more succinct, the exposure shift defines a
counterfactual scenario. For this application, the threshold $c$ is evaluated
at equally spaced points starting with 15 $\mathrm{\upmu g/m^{3}}$ moving down
to 6 $\mathrm{\upmu g/m^{3}}$. We expect that at 15 $\mathrm{\upmu g/m^{3}}$
there will be little to no reduction in deaths since $>99\%$ of observations
fall below that range. We can contemplate the proposed NAAQS levels through
$d$ assuming that full compliance to the new regulation holds for incompliant
ZIP-codes. The exposure shift should otherwise not affect already compliant
ZIP-codes. Second, we onsider _percent reduction shifts_. This scenario
assumes that all ZIP-code years reduce their pollution levels proportionally
from their observed value. More precisely, the shift is defined as
$\tilde{A}=A(1-c)$. We considered a range of percent reduction shifts between
$c\in(0,50)\%$. We can interpret these shifts in terms of the NAAQS by mapping
each percent reduction to a compliance percentile. For instance, Figure 4
shows that, under a 30% overall reduction in historical values, approximately
82% would comply with a NAAQS of 9 $\mathrm{\upmu g/m^{3}}$.

Implementation We implement tresnet using varying-coefficient splines as in
Section 5. We select a NN architecture using the out-of-sample prediction
error from a 20/80% test-train split to choose the number of hidden layers
(1-3 layers) and hidden dimensions (16, 64, 256). We found no evidence of
overfitting in the selected models. To account for uncertainty in our
estimations, we train the model on 100 bootstrap samples, each with random
initializations, thereby obtaining an approximate posterior distribution. Deep
learning ensembles have been previously shown to approximate estimation
uncertainty well in deep learning tasks (Izmailov et al., 2021).

Figure 5: Estimated SRF of the total deaths (%) for different cutoffs.

Results Figure 1 in the introduction presents the effects of shifting the
$\text{PM}_{2.5}$ distribution at various cutoffs on the expected reduction to
deaths. The slope is steeper at stricter/lower cutoffs, likely because lower
cutoffs affect a larger fraction of the observed population and reduce the
overall $\text{PM}_{2.5}$. For instance, figure Figure 1 shows that had no
ZIP-code years exceeded 12 $\mathrm{\upmu g/m^{3}}$, the observed death counts
would have decreased by around 1%. If the cutoff is lowered to 9
$\mathrm{\upmu g/m^{3}}$, then deaths could have fallen by around 4%. The
slope becomes increasingly steeper as the $\text{PM}_{2.5}$ threshold is
reduced, suggesting the increasing benefits of lowering the standard
concentration level. Another way to interpret this result is to say that there
is a greater gain to reducing mortality caused by $\text{PM}_{2.5}$ from
lowering the concentration level from 10 to 8 $\mathrm{\upmu g/m^{3}}$ than
there is from lowering it from 12 to 10 $\mathrm{\upmu g/m^{3}}$.

The results of the percent-reduction shift are presented in Figure 5. The
decrease in deaths is approximately linear with respect to the percent
decrease in $\text{PM}_{2.5}$. As such, the SRF shows an approximate 0.5%
decrease in deaths resulting from a 10% decrease in $\text{PM}_{2.5}$. This
result is consistent with previous causal estimates of the marginal effect of
$\text{PM}_{2.5}$ exposure on elder mortality (Wu et al., 2020). Percent
reduction offers a complementary view to the cutoff shift response function.

##  7 Discussion and Limitations

We have made a significant stride in addressing the pressing public health
question regarding the potential health benefits of lowering the NAAQS in the
United States. In response to this question, we introduce the first causal
inference method to utilize neural networks for estimating SRFs. Furthermore,
we have extended this method to handle count data, which is crucial for our
application in addition to other public health and epidemiology contexts. We
acknowledge some limitations to our methodology. First, our uncertainty
assessment of the SRF relies on the bootstrap and ensembling of multiple
random seeds. While these methods are used often in practice, future research
could explore the integration of tresnet with Bayesian methods to enhance
uncertainty quantification. Second, our application of the methodology focuses
on exposure shifts representing complementary viewpoints to the possible
effects of the proposed EPA rules on the NAAQS. However, it does not determine
the most probable exposure shift resulting from the new rule’s implementation,
based on historical responses to changes in the NAAQS. Subsequent
investigations should more carefully consider this aspect of the analysis. The
assessment of annual average $\text{PM}_{2.5}$ levels at the ZIP-code level is
based on predictions rather than on actual observable values, introducing
potential attenuation bias stemming from measurement error. Nonetheless,
previous studies on measurement error involving clustered air pollution
exposures have demonstrated that such attenuation tends to pull the causal
effect towards a null result (Josey et al., 2022; Wei et al., 2022). It is
essential to recognize that the SRF framework places additional considerations
on the analyst designing the exposure shift. This newfound responsibility can
be seen as both a disadvantage and an advantage. However, it highlights the
need for an explicit and meticulous statement of the assumptions underlying
the considered exposure shifts in order to mitigate the potential misuse of
SRF estimation techniques.

## References

  * Bahadori et al. (2022) Taha Bahadori, Eric Tchetgen Tchetgen, and David Heckerman.  End-to-end balancing for causal continuous treatment-effect estimation.  In _International Conference on Machine Learning_ , pp. 1313–1326, 2022.
  * Bang & Robins (2005) Heejung Bang and James M Robins.  Doubly robust estimation in missing data and causal inference models.  _Biometrics_ , 61(4):962–973, 2005.
  * Bica et al. (2020) Ioana Bica, James Jordon, and Mihaela van der Schaar.  Estimating the effects of continuous-valued interventions using generative adversarial networks.  _Advances in Neural Information Processing Systems_ , 33:16434–16445, 2020.
  * Bickel et al. (1993) Peter J Bickel, Chris AJ Klaassen, Peter J Bickel, Ya’acov Ritov, J Klaassen, Jon A Wellner, and YA’Acov Ritov.  _Efficient and adaptive estimation for semiparametric models_ , volume 4.  Springer, 1993.
  * Chiang et al. (2001) Chin-Tsang Chiang, John A Rice, and Colin O Wu.  Smoothing spline estimation for varying coefficient models with repeatedly measured dependent variables.  _Journal of the American Statistical Association_ , 96(454):605–619, 2001.
  * Di et al. (2019) Qian Di, Heresh Amini, Liuhua Shi, Itai Kloog, Rachel Silvern, James Kelly, M Benjamin Sabath, Christine Choirat, Petros Koutrakis, Alexei Lyapustin, et al.  An ensemble-based model of pm2. 5 concentration across the contiguous united states with high spatiotemporal resolution.  _Environment international_ , 130:104909, 2019.
  * Díaz & Hejazi (2020) Iván Díaz and Nima S Hejazi.  Causal mediation analysis for stochastic interventions.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 82(3), 2020.
  * Díaz et al. (2021) Iván Díaz, Nicholas Williams, Katherine L Hoffman, and Edward J Schenck.  Nonparametric causal effects based on longitudinal modified treatment policies.  _Journal of the American Statistical Association_ , pp. 1–16, 2021\.
  * Farrell et al. (2021) Max H Farrell, Tengyuan Liang, and Sanjog Misra.  Deep neural networks for estimation and inference.  _Econometrica_ , 89(1):181–213, 2021.
  * Hill (2011) Jennifer L Hill.  Bayesian nonparametric modeling for causal inference.  _Journal of Computational and Graphical Statistics_ , 20(1):217–240, 2011.
  * Imbens & Rubin (2015) Guido W. Imbens and Donald B. Rubin.  _Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction_.  Cambridge University Press, 2015.
  * Izmailov et al. (2021) Pavel Izmailov, Sharad Vikram, Matthew D Hoffman, and Andrew Gordon Gordon Wilson.  What are bayesian neural network posteriors really like?  In _International conference on machine learning_ , pp. 4629–4640, 2021.
  * Josey et al. (2022) Kevin P Josey, Priyanka DeSouza, Xiao Wu, Danielle Braun, and Rachel Nethery.  Estimating a causal exposure response function with a continuous error-prone exposure: A study of fine particulate matter and all-cause mortality.  _Journal of Agricultural, Biological and Environmental Statistics_ , pp. 1–22, 2022.
  * Josey et al. (2023) Kevin P Josey, Scott W Delaney, Xiao Wu, Rachel C Nethery, Priyanka DeSouza, Danielle Braun, and Francesca Dominici.  Air pollution and mortality at the intersection of race and social class.  _New England Journal of Medicine_ , 388(15), 2023.
  * Kennedy (2016) Edward H Kennedy.  Semiparametric theory and empirical processes in causal inference.  _Statistical causal inferences and their applications in public health research_ , pp. 141–167, 2016.
  * Kennedy (2022) Edward H Kennedy.  Semiparametric doubly robust targeted double machine learning: a review.  _arXiv preprint arXiv:2203.06469_ , 2022.
  * McCullagh (2019) Peter McCullagh.  _Generalized linear models_.  Routledge, 2019.
  * Muñoz & Van Der Laan (2012) Iván Díaz Muñoz and Mark Van Der Laan.  Population intervention causal effects based on stochastic interventions.  _Biometrics_ , 68(2):541–549, 2012.
  * Newman (2008) David Newman.  Bag of words data set.  _UCI Machine Learning Respository_ , 289, 2008.
  * Nie et al. (2021) Lizhen Nie, Mao Ye, Dan Nicolae, et al.  Vcnet and functional targeted regularization for learning causal effects of continuous treatments.  In _International Conference on Learning Representations_ , 2021.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32, 2019.
  * Robins (2000) James M Robins.  Robust estimation in sequentially ignorable missing data and causal inference models.  In _Proceedings of the American Statistical Association_ , volume 1999, 2000.
  * Robins et al. (2000) James M Robins, Andrea Rotnitzky, and Mark van der Laan.  On profile likelihood: comment.  _Journal of the American Statistical Association_ , 95(450):477–482, 2000.
  * Schwab et al. (2020) Patrick Schwab, Lorenz Linhardt, Stefan Bauer, Joachim M Buhmann, and Walter Karlen.  Learning counterfactual representations for estimating individual dose-response curves.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pp. 5612–5619, 2020.
  * Shi et al. (2019) Claudia Shi, David Blei, and Victor Veitch.  Adapting neural networks for the estimation of treatment effects.  _Advances in neural information processing systems_ , 32, 2019.
  * Smith et al. (2023) Matthew J Smith, Rachael V Phillips, Miguel Angel Luque-Fernandez, and Camille Maringe.  Application of targeted maximum likelihood estimation in public health and epidemiological studies: a systematic review.  _Annals of Epidemiology_ , 2023.
  * Tsiatis (2006) Anastasios A Tsiatis.  Semiparametric theory and missing data.  2006\.
  * Van der Laan et al. (2011) Mark J Van der Laan, Sherri Rose, et al.  _Targeted learning: causal inference for observational and experimental data_ , volume 4.  Springer, 2011.
  * Wainwright (2019) Martin J Wainwright.  _High-dimensional statistics: A non-asymptotic viewpoint_ , volume 48.  Cambridge university press, 2019.
  * Wei et al. (2022) Yaguang Wei, Xinye Qiu, Mahdieh Danesh Yazdi, Alexandra Shtein, Liuhua Shi, Jiabei Yang, Adjani A Peralta, Brent A Coull, and Joel D Schwartz.  The impact of exposure measurement error on the estimated concentration–response relationship between long-term exposure to pm 2.5 and mortality.  _Environmental Health Perspectives_ , 130(7):077006, 2022.
  * Wu et al. (2020) X Wu, D Braun, J Schwartz, MA Kioumourtzoglou, and F Dominici.  Evaluating the impact of long-term exposure to fine particulate matter on mortality among the elderly.  _Science advances_ , 6(29):eaba5692, 2020.
  * Yoon et al. (2018) Jinsung Yoon, James Jordon, and Mihaela Van Der Schaar.  Ganite: Estimation of individualized treatment effects using generative adversarial nets.  In _International conference on learning representations_ , 2018.

##  Appendix A Background on the EIF

When we fit a statistical model to a dataset, each data point contributes to
the estimated parameters of the model. The efficient influence function (EIF)
effectively measures how sensitive an estimate is to the inclusion or
exclusion of the individual data points. In other words, it tells us how much
a single data point can influence the estimate of a causal effect. The term
“efficient” refers to the property that, over the aggregate, the efficient
influence function evaluated over the observed data provides the uniformly
best statistically-efficient estimator of the associated causal effect under
the nonparametric identifying assumptions that were made explicit in the main
body of the manuscript (Kennedy, 2016). This property holds even when using
machine learning methods (like neural networks) to estimate the nuisance
parameters – in our case these nuisance parameters refer to the outcome
regression model and the generalized propensity score model (Muñoz & Van Der
Laan, 2012). An important distinction to make in this definition is that
efficiency is an asymptotic property, meaning that it only holds as the sample
size goes to infinity.

Another way to describe the efficient influence function is as the canonical
gradient (Van der Laan et al., 2011) of the targeted parameter. Specifically,
let $\\{{\mathbb{P}}_{t}:t\in\mathbb{R}\\}$ be a smooth parametric submodel
such that ${\mathbb{P}}_{0}={\mathbb{P}}$ and denote its score function as
$s_{h}({\bm{o}})=\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\big{|}_{t=h}\log\tilde{p}_{t}({\bm{o}})$
where we use
$\tilde{p}_{t}({\bm{o}})=p_{t}(y|{\bm{x}},a)\tilde{p}_{t}(a|{\bm{x}})p_{t}({\bm{x}})$
and $p_{t}({\bm{o}})=p_{t}(y|{\bm{x}},a){p}_{t}(a|{\bm{x}})p_{t}({\bm{x}})$ to
denote the probability density function of ${\mathbb{P}}_{t}$. We also define
$w_{t}({\bm{x}},a)=\tilde{p}_{t}(a|{\bm{x}})/p_{t}(a|{\bm{x}})$. The SRF
estimand at each member of the submodel is defined as

$\displaystyle\psi({\mathbb{P}}_{t})=\int\int\int
y\tilde{p}_{t}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}.$

For (9) to be the efficient influence function in a non-parametric model, one
must show that it is _pathwise differentiable_ , that is,

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\Big{|}_{t=0}\psi({\mathbb{P}}_{t})=\int\varphi({\bm{o}};\psi,\mu,{\bm{w}})s_{0}({\bm{o}})\tilde{p}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}.$
(8)

This proof appears in the next section.

##  Appendix B Technical Proofs

Before proving the main theorems in the paper, we show the following result
regarding the identification of the ${\bm{\psi}}$.

Suppose Figure 2 holds. Then $\mu({\bm{X}},A)=\mathbb{E}[Y\mid{\bm{X}},A]$.
The right-hand side does not involve potential outcomes. As a corollary,
${\bm{\psi}}=\mathbb{E}[\mu({\bm{X}},\tilde{A})]$ can be identified.

###### Proof.

Since the treatment and potential outcomes are independent conditional on
${\bm{X}}$ (unconfoundedness). By properties of the conditional expectation,
we have

$\displaystyle\mu({\bm{x}},a)$
$\displaystyle=\mathbb{E}[Y^{a}\mid{\bm{X}}={\bm{x}}]$
$\displaystyle=^{(a)}\mathbb{E}[Y^{a}\mid{\bm{X}}={\bm{x}},A=a]$
$\displaystyle=\mathbb{E}[Y^{A}\mid{\bm{X}}={\bm{x}},A=a]=^{(b)}\mathbb{E}[Y\mid{\bm{X}}={\bm{x}},A=a].$

The identity (a) uses unconfoundedness (Figure 2); (b) uses that $Y^{A}=Y$
(consistency). ∎

We now prove the three main theorems. The proof strategy follows the same
pattern as Nie et al. (2021), with modifications to accommodate the SRF
estimand and is mindful of the exponential family loss, which complicates the
analytical form of the fluctuation term.

[] Suppose Assumptions 2 and 2 hold. Then the EIF of ${\bm{\psi}}$ is given by

$\displaystyle\varphi({\bm{O}};{\bm{\psi}},\mu,{\bm{w}})={\bm{w}}({\bm{X}},A)\left(Y-\mu({\bm{X}},A)\right)+\mu({\bm{X}},\tilde{A})-{\bm{\psi}}.$
(9)

Furthermore, let $\hat{\mu}=\hat{\mu}({\mathbb{P}}_{n})$ and
$\hat{{\bm{w}}}=\hat{{\bm{w}}}({\mathbb{P}}_{n})$ be estimators such that
$\lVert\hat{\mu}-\mu\rVert_{\infty}=O_{p}(r_{1}(n))$ and
$\lVert\hat{{\bm{w}}}-{\bm{w}}\rVert_{\infty}=O_{p}(r_{2}(n))$. Then
$\lVert{\mathbb{P}}\varphi({\bm{\psi}},\hat{\mu},\hat{{\bm{w}}})\rVert=O_{p}(r_{1}(n)r_{2}(n))$.
One of the conditions $\hat{\mu}=\mu$ or $\hat{{\bm{w}}}={\bm{w}}$ suffice for
${\mathbb{P}}\varphi({\bm{\psi}},\hat{\mu},\hat{{\bm{w}}})={\bm{0}}$ to hold.

### Step 1: showing that $\varphi$ is the EIF

###### Proof.

Recall the notation from Appendix A. To show that
$\varphi({\bm{o}};{\bm{\psi}},\mu,{\bm{w}})$ in Equation 9 is the EIF, we must
show that

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\Big{|}_{t=0}\psi({\mathbb{P}}_{t})=\int\varphi({\bm{o}};\psi,\mu,{\bm{w}})s_{0}({\bm{o}})\tilde{p}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
(10)

Due to the logarithmic transformation defining the score equation, we can
factorize $s_{t}({\bm{o}})$ into three parts with
$s_{t}({\bm{o}})=s_{1,t}(y|a,{\bm{x}})+s_{2,t}(a|{\bm{x}})+s_{3,t}({\bm{x}})$.

Starting on the left hand side of (10), under standard regularity conditions
we have:

$\displaystyle\int
ys_{t}({\bm{o}})\tilde{p}_{t}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
$\displaystyle=\int
y\frac{\mathop{}\\!\mathrm{d}{\tilde{p}_{t}({\bm{o}})}/\mathop{}\\!\mathrm{d}{t}}{\tilde{p}_{t}({\bm{o}})}p_{t}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
(11) $\displaystyle=\int
y\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}{t}}{\tilde{p}_{t}({\bm{o}})}\mathop{}\\!\mathrm{d}{{\bm{o}}}$
$\displaystyle=\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}{t}}\int
y\tilde{p}_{t}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}},$

which allows us to partition the pathwise derivative (i.e. the left-hand side
of (10)) into:

$\displaystyle\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}{t}}\psi({\mathbb{P}}_{t})$
$\displaystyle=\int
ys_{t}(\mathbf{o})\tilde{p}_{t}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
(12) $\displaystyle=\int
y(s_{1,t}(y|a,{\bm{x}})+s_{2,t}(a|{\bm{x}})+s_{3,t}({\bm{x}}))\tilde{p}_{t}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
$\displaystyle=\int
ys_{1,t}(y|a,{\bm{x}})\tilde{p}_{t}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}+\int
ys_{2,t}(a|{\bm{x}})\tilde{p}_{t}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}+\int
ys_{3,t}({\bm{x}})\tilde{p}_{t}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}.$

For the next step of the proof, we will recursively use the two following
identities. For any arbitrary function $g(\cdot)$, and for any two subsets of
measurements ${\bm{o}}_{1}$ and ${\bm{o}}_{2}$ (e.g.
${\bm{o}}_{1}=(a,{\bm{x}})$ and ${\bm{o}}_{2}=y$), we have

$\displaystyle\int
g({\bm{o}}_{1})\left\\{{\bm{o}}_{2}-\int{\bm{o}}_{2}\tilde{p}_{t}({\bm{o}}_{2}|{\bm{o}}_{1})\mathop{}\\!\mathrm{d}{{\bm{o}}_{2}}\right\\}\tilde{p}_{t}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
$\displaystyle=^{(a)}0$ (13) $\displaystyle\int
g({\bm{o}}_{1})s_{t}({\bm{o}}_{2}|{\bm{o}}_{1})\tilde{p}_{t}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
$\displaystyle=^{(b)}0$

Defining $\tilde{p}_{0}(\cdot)=\tilde{p}(\cdot)$ and $p_{0}(\cdot)=p(\cdot)$,
we can then find the first term in (12) is equal to

$\displaystyle\int
ys_{1,t}(y|a,{\bm{x}})\mathop{}\\!\mathrm{d}\tilde{p}_{t}({\bm{o}})\mathop{}\\!\mathrm{d}{\bm{o}}\Big{|}_{t=0}$
$\displaystyle=\int
ys_{1,0}(y|a,{\bm{x}})\tilde{p}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
$\displaystyle=^{(b)}\int\left\\{y-\mu({\bm{x}},a)\right\\}s_{1,0}(y|a,{\bm{x}})p(y|a,{\bm{x}})\tilde{p}(a|{\bm{x}})p({\bm{x}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
$\displaystyle=\int
w_{t}({\bm{x}},a)\left\\{y-\mu({\bm{x}},a)\right\\}s_{1,0}(y|a,{\bm{x}})p(y|a,{\bm{x}})p(a|{\bm{x}})p({\bm{x}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
$\displaystyle=^{(a)}\int
w_{t}({\bm{x}},a)\left\\{y-\mu({\bm{x}},a)\right\\}\left\\{s_{1,0}(y|a,{\bm{x}})+s_{2,0}(a|{\bm{x}})\right\\}p({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
$\displaystyle=^{(a)}\int
w_{t}({\bm{x}},a)\left\\{y-\mu({\bm{x}},a)\right\\}\left\\{s_{1,0}(y|a,{\bm{x}})+s_{2,0}(a|{\bm{x}})+s_{3,0}({\bm{x}})\right\\}p({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
$\displaystyle=\int
w_{t}({\bm{x}},a)\left\\{y-\mu({\bm{x}},a)\right\\}s_{0}({\bm{o}})p({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}.$

For the second term, we have

$\displaystyle\int
ys_{2,t}(a|{\bm{x}})\mathop{}\\!\mathrm{d}\tilde{p}_{t}({\bm{o}})\mathop{}\\!\mathrm{d}{\bm{o}}\Big{|}_{t=0}$
$\displaystyle=\int\int\int
ys_{2,0}(a|{\bm{x}})p(y|a,{\bm{x}})\tilde{p}(a|{\bm{x}})p({\bm{x}})\mathop{}\\!\mathrm{d}{y}\mathop{}\\!\mathrm{d}{a}\mathop{}\\!\mathrm{d}{{\bm{x}}}$
$\displaystyle=\int\int\mu({\bm{x}},a)s_{2,0}(a|{\bm{x}})\tilde{p}(a|{\bm{x}})p({\bm{x}})\mathop{}\\!\mathrm{d}{a}\mathop{}\\!\mathrm{d}{{\bm{x}}}$
$\displaystyle=^{(b)}\int\int\left\\{\mu({\bm{x}},a)-\int\mu({\bm{x}},a)\tilde{p}(a|{\bm{x}})\mathop{}\\!\mathrm{d}{a}\right\\}s_{2,0}(a|{\bm{x}})\tilde{p}(a|{\bm{x}})p({\bm{x}})\mathop{}\\!\mathrm{d}{a}\mathop{}\\!\mathrm{d}{{\bm{x}}}$
$\displaystyle=^{(b)}\int\left\\{\mu({\bm{x}},a)-\int\mu({\bm{x}},a)\tilde{p}(a|{\bm{x}})\mathop{}\\!\mathrm{d}{a}\right\\}\left\\{s_{1,0}(y|a,{\bm{x}})+s_{2,0}(a|{\bm{x}})\right\\}\tilde{p}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
$\displaystyle=^{(a)}\int\left\\{\mu({\bm{x}},a)-\int\mu({\bm{x}},a)\tilde{p}(a|{\bm{x}})\mathop{}\\!\mathrm{d}{a}\right\\}\left\\{s_{1,0}(y|a,{\bm{x}})+s_{2,0}(a|{\bm{x}})+s_{3,0}({\bm{x}})\right\\}\tilde{p}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
$\displaystyle=\int\left\\{\mu({\bm{x}},a)-\int\mu({\bm{x}},a)\tilde{p}(a|{\bm{x}})\mathop{}\\!\mathrm{d}{a}\right\\}s_{0}({\bm{o}})\tilde{p}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$

where we again center the integrand and apply the identities of (13). Finally,
for the third term in (12) we have

$\displaystyle\int
ys_{3,t}(\mathbf{x})\mathop{}\\!\mathrm{d}\tilde{p}_{t}({\bm{o}})\mathop{}\\!\mathrm{d}{\bm{o}}\Big{|}_{t=0}$
$\displaystyle=\int\int\int
ys_{3,0}({\bm{x}})p(y|a,{\bm{x}})\tilde{p}(a|{\bm{x}})p({\bm{x}})\mathop{}\\!\mathrm{d}{y}\mathop{}\\!\mathrm{d}{a}\mathop{}\\!\mathrm{d}{x}$
$\displaystyle=\int\left\\{\int\mu({\bm{x}},a)\tilde{p}(a|{\bm{x}})\mathop{}\\!\mathrm{d}{a}\right\\}s_{3,0}({\bm{x}})p({\bm{x}})\mathop{}\\!\mathrm{d}{x}$
$\displaystyle=^{(b)}\int\left\\{\int\mu({\bm{x}},a)\tilde{p}(a|{\bm{x}})\mathop{}\\!\mathrm{d}{a}-\psi\right\\}s_{3,0}({\bm{x}})p({\bm{x}})\mathop{}\\!\mathrm{d}{x}$
$\displaystyle=^{(b)}\int\left\\{\int\mu({\bm{x}},a)\tilde{p}(a|{\bm{x}})\mathop{}\\!\mathrm{d}{a}-\psi\right\\}\left\\{s_{1,0}(y|a,{\bm{x}})+s_{2,0}(a|{\bm{x}})+s_{3,0}({\bm{x}})\right\\}\tilde{p}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$
$\displaystyle=\int\left\\{\int\mu({\bm{x}},a)\tilde{p}(a|{\bm{x}})\mathop{}\\!\mathrm{d}{a}-\psi\right\\}s_{0}({\bm{o}})\tilde{p}({\bm{o}})\mathop{}\\!\mathrm{d}{{\bm{o}}}$

Combining these three terms above proves the condition in (10) holds, thus
completing the proof. ∎

### Step 2: showing the estimating equation and double robustness

###### Proof.

We can break down ${\mathbb{P}}\varphi$ in two terms,

${\mathbb{P}}\varphi({\bm{\psi}},\hat{\mu},\hat{{\bm{w}}})=\mathbb{E}[\hat{{\bm{w}}}({\bm{X}},A)(Y-\hat{\mu}({\bm{X}},A))]+\mathbb{E}[\hat{\mu}(X,\tilde{A})-{\bm{\psi}}]$

We can rewrite the first term as

$\mathbb{E}[\hat{{\bm{w}}}({\bm{X}},A)(Y-\hat{\mu}({\bm{X}},A))]=\mathbb{E}[\hat{{\bm{w}}}({\bm{X}},A)(\mu({\bm{X}},A)-\hat{\mu}({\bm{X}},A))]$

For the second term, we can write

$\displaystyle\mathbb{E}[\hat{\mu}({\bm{X}},\tilde{A})-{\bm{\psi}}]$
$\displaystyle=\mathbb{E}[\hat{\mu}-\mu({\bm{X}},\tilde{A})]$ (14)
$\displaystyle=\mathbb{E}[{\bm{w}}({\bm{X}},A)(\hat{\mu}({\bm{X}},\tilde{A})-\mu({\bm{X}},\tilde{A}))]$

Combining them, we have that

${\mathbb{P}}\varphi({\bm{\psi}},\hat{\mu},\hat{{\bm{w}}})=\mathbb{E}[(\hat{{\bm{w}}}({\bm{X}},A)-{\bm{w}}({\bm{X}},A))(\mu({\bm{X}},A)-\hat{\mu}({\bm{X}},A))].$

From this expression, the conclusion of the theorem is evident.
${\mathbb{P}}\varphi({\bm{\psi}},\hat{\mu},\hat{{\bm{w}}})={\bm{0}}$ if either
$\hat{\mu}=\mu$ or $\hat{{\bm{w}}}={\bm{w}}$. Further,

$\lVert\mathbb{E}[(\hat{{\bm{w}}}({\bm{X}},A)-{\bm{w}}({\bm{X}},A))(\mu({\bm{X}},A)-\hat{\mu}({\bm{X}},A))]\rVert\leq\lVert\hat{{\bm{w}}}-{\bm{w}}\rVert_{\infty}\lVert\hat{\mu}-\mu\rVert_{\infty}=O_{p}(r_{1}(n)r_{2}(n)),$

completing the proof. ∎

We can now prove the first theorem. In what follows we will adopt the notation
${\mathbb{P}}_{n}$ to denote the empirical distribution. For any function
$f({\bm{O}})$, we have that ${\mathbb{P}}_{n}f=\frac{1}{n}f(O_{i})$ and
${\mathbb{P}}f=\mathbb{E}[f({\bm{O}})]$.

See 3

###### Proof.

A key property of the exponential family of distributions and the associated
link function $g$ that we will use are
$\Lambda^{\prime}(\eta)=\frac{d}{d\eta}\Lambda(\eta)=g^{-1}(E[Y|\eta])$
(McCullagh, 2019). Also note that
$\frac{d}{d{\bm{\epsilon}}}g(\tilde{\mu}^{\text{NN}}({\bm{X}},A))={\bm{w}}^{\text{NN}}({\bm{X}},A)$
for all $\mu^{\text{NN}},{\bm{w}}^{\text{NN}}{\bm{\epsilon}}$. These two
observations and the chain rule give us

$\displaystyle{\bm{0}}$
$\displaystyle=\frac{d}{d{\bm{\epsilon}}}{\mathcal{R}}^{\text{tr}}(\hat{\mu},\hat{{\bm{w}}},\hat{{\bm{\epsilon}}})$
(15)
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}{\bm{\epsilon}}}\Big{|}_{{\bm{\epsilon}}=\hat{{\bm{\epsilon}}}}\left\\{\Lambda(g(\tilde{\mu}({\bm{X}},A)))-Yg(\tilde{\mu}({\bm{X}},A))\right\\}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left\\{g^{-1}(g(\tilde{\mu}({\bm{X}},A)))\hat{{\bm{w}}}({\bm{X}},A)-Y\hat{{\bm{w}}}({\bm{X}},A)\right\\}=\frac{1}{n}\sum_{i=1}^{n}\hat{{\bm{w}}}({\bm{X}},A)(\tilde{\mu}({\bm{X}},A)-Y).$

The fact that
$\hat{{\bm{\psi}}}^{\text{tr}}=\frac{1}{n}\sum_{i=1}^{n}\tilde{\mu}({\bm{X}}_{i},\tilde{A}_{i})$
satisfies the empirical estimating equation follows trivially from the fact
that
${\mathbb{P}}_{n}\varphi(\hat{{\bm{\psi}}}^{\text{tr}},\tilde{\mu},\hat{{\bm{w}}})=\frac{1}{n}\sum_{i=1}^{n}\hat{{\bm{w}}}({\bm{X}},A)(Y-\tilde{\mu}({\bm{X}},A))+\frac{1}{n}\sum_{i=1}^{n}\tilde{\mu}({\bm{X}}_{i},\tilde{A}_{i})-\hat{{\bm{\psi}}}^{\text{tr}}$.
The first term is zero because of the above results while the last two terms
cancel each other by definition. ∎

See 3

### Step 1: showing consistency of the outcome and density ratio models

This proof closely adapts Nie et al. (2021) for the ERF case. We will use the
notation
${\mathcal{R}}_{\mathbb{P}}(\mu^{\text{NN}},{\bm{w}}^{\text{NN}})={\mathbb{P}}{\mathcal{L}}(\mu^{\text{NN}},{\bm{w}}^{\text{NN}})$
for the population risk, where ${\mathcal{L}}^{\text{NN}}$ is the loss
function. Denote $\mu^{*},{\bm{w}}^{*}$ as the population risk minimizers. We
assume $\mu^{*}\in{\mathcal{M}},{\bm{w}}^{*}\in{\mathcal{W}}$. The proof
strategy is fairly standard and requires that the loss function be Lipschitz
to ensure that the Rademacher complexity and boundedness assumptions extend
the loss terms. The Lipschitz condition can often be relaxed with direct
assumptions on the terms ${\mathcal{L}}(\mu^{\text{NN}},{\bm{w}}^{\text{NN}})$
to have vanishing Rademacher complexity and boundedness (Wainwright, 2019).

###### Proof.

We first show that the risk of the regularized parameters is not too different
from the minimum population risk. Specifically, we show that

${\mathcal{R}}_{\mathbb{P}}^{\text{NN}}(\hat{\mu},\hat{{\bm{w}}})-{\mathcal{R}}_{\mathbb{P}}^{\text{NN}}(\mu^{*},{\bm{w}}^{*})=o(1)+O_{p}(n^{-1/2}).$
(16)

To prove this fact, we first note that

$\displaystyle 0$
$\displaystyle\leq{\mathcal{R}}_{\mathbb{P}}^{\text{NN}}(\hat{\mu},\hat{{\bm{w}}})-{\mathcal{R}}_{\mathbb{P}}^{\text{NN}}(\mu^{*},{\bm{w}}^{*})$
(17)
$\displaystyle\leq{\mathcal{R}}(\hat{\mu},\hat{{\bm{w}}})-{\mathcal{R}}^{\text{NN}}(\mu^{*},{\bm{w}}^{*}))+({\mathbb{P}}-{\mathbb{P}}_{n}){\mathcal{L}}^{\text{NN}}(\hat{\mu},\hat{{\bm{w}}})+({\mathbb{P}}_{n}-{\mathbb{P}}){\mathcal{L}}^{\text{NN}}(\mu^{*},{\bm{w}}^{*}).$

The second and third terms are empirical processes. We now use the regularity
assumptions on the Rademacher complexity. The order of the Rademacher
complexity is preserved under Lipschitz transforms. Hence, under the
assumption that the loss is Lipschitz, it follows that the class
$\\{{\mathcal{L}}^{\text{NN}}(\mu^{\text{NN}},{\bm{w}}^{\text{NN}})\colon\mu^{\text{NN}}\in{\mathcal{M}},{\bm{w}}^{\text{NN}}\in{\mathcal{W}}\\}$
has a Rademacher complexity of order $O(n^{-1/2})$. Uniform boundedness is
also preserved under Lispschitz transformations. Together, the vanishing
Rademacher complexity and uniform boundedness imply the uniform law of large
numbers, which in turn implies the convergence of the empirical process
(Wainwright, 2019). Hence the last two terms are $O_{p}(n^{-1/2})$. We now
bound the first term.

$\displaystyle{\mathcal{R}}^{\text{NN}}(\hat{\mu},\hat{{\bm{w}}})-{\mathcal{R}}^{\text{NN}}(\mu^{*},{\bm{w}}^{*})$
(18)
$\displaystyle\quad=({\mathcal{R}}^{\text{NN}}+\beta_{n}{\mathcal{R}}^{\text{tr}})(\hat{\mu},\hat{{\bm{w}}},\hat{{\bm{\epsilon}}})-({\mathcal{R}}^{\text{NN}}+\beta_{n}{\mathcal{R}}^{\text{tr}})(\mu^{*},{\bm{w}}^{*},{\bm{0}})+\beta_{n}({\mathcal{R}}^{\text{tr}}(\mu^{*},{\bm{w}}^{*},{\bm{0}})-{\mathcal{R}}^{\text{tr}}(\hat{\mu},\hat{{\bm{w}}},\hat{{\bm{\epsilon}}}))$
$\displaystyle\quad\leq^{(a)}\beta_{n}({\mathcal{R}}^{\text{tr}}(\mu^{*},{\bm{w}}^{*},{\bm{0}})-{\mathcal{R}}^{\text{tr}}(\hat{\mu},\hat{{\bm{w}}},\hat{{\bm{\epsilon}}}))$
$\displaystyle\quad=^{(b)}\beta_{n}({\mathbb{P}}_{n}{\mathcal{L}}_{\mu}(\mu^{*})+O(1)).$
$\displaystyle\quad=\beta_{n}(({\mathbb{P}}_{n}-{\mathbb{P}}){\mathcal{L}}_{\mu}(\mu^{*})+{\mathbb{P}}{\mathcal{L}}_{\mu}(\mu^{*})+O(1))$
$\displaystyle\quad=^{(c)}\beta_{n}(O_{p}(n^{-1/2}))+O(1))=o_{p}(1)$

Inequality (a) is due to $(\hat{\mu},\hat{{\bm{w}}},\hat{{\bm{\epsilon}}})$
being a minimizer for the regularized risk. Inequality (b) is the result of
${\mathcal{R}}^{\text{tr}}(\hat{\mu},\hat{{\bm{w}}},\hat{{\bm{\epsilon}}})$
being bounded since the exponential family of distributions is log-concave and
${\mathcal{R}}^{\text{tr}}(\mu^{*},{\bm{w}}^{*},{\bm{0}})={\mathbb{P}}_{n}{\mathcal{L}}_{\mu}(g(\mu^{*}))$.
(c) uses the uniform law of large numbers from the Rademacher complexity and
the uniform boundedness, and the fact that ${\mathcal{L}}$ is Lipschitz.
Combining Equation 17 and Equation 18, we get
${\mathcal{R}}_{\mathbb{P}}^{\text{NN}}(\hat{\mu},\hat{{\bm{w}}})-{\mathcal{R}}_{\mathbb{P}}^{\text{NN}}(\mu^{*},{\bm{w}}^{*})=o_{p}(1)$.

The result now follows from observing that the population risk has a unique
minimizer up to the reparameterization of the network weights. Hence, by
regularity conditions, $\lVert\hat{\mu}-\mu\rVert=o_{p}(1)$ and
$\lVert\hat{{\bm{w}}}-{\bm{w}}\rVert=o_{p}(1)$. ∎

### Step 2: Proving convergence and efficiency of
$\hat{{\bm{\psi}}}^{\text{tr}}$

###### Proof.

Direct computation gives

$\displaystyle\lVert\hat{{\bm{\psi}}}^{\text{tr}}-{\bm{\psi}}\rVert$
$\displaystyle=\lVert\textstyle{\frac{1}{n}\sum_{i=1}^{n}\tilde{\mu}({\bm{X}}_{i},\tilde{A}_{i})}-{\bm{\psi}}\rVert$
(19)
$\displaystyle=^{(a)}\lVert\textstyle{\frac{1}{n}\sum_{i=1}^{n}\\{\tilde{\mu}({\bm{X}}_{i},\tilde{A}_{i})+\hat{{\bm{w}}}({\bm{X}}_{i},A_{i})(Y_{i}-\tilde{\mu}({\bm{X}}_{i},A_{i}))\\}}-{\bm{\psi}}\rVert$
$\displaystyle=^{(b)}\lVert\mathbb{E}[\hat{{\bm{w}}}({\bm{X}},A)(Y-\tilde{\mu}({\bm{X}},A))+\tilde{\mu}({\bm{X}},\tilde{A})]-{\bm{\psi}}\rVert+O_{p}(n^{-1/2}),$
$\displaystyle=^{(c)}\lVert\mathbb{E}[\hat{{\bm{w}}}({\bm{X}},A)(\mu({\bm{X}},A)-\tilde{\mu}({\bm{X}},A))+\tilde{\mu}({\bm{X}},\tilde{A})]-{\bm{\psi}}\rVert+O_{p}(n^{-1/2}),$
$\displaystyle=^{(d)}\lVert\mathbb{E}[(\hat{{\bm{w}}}({\bm{X}},A)-{\bm{w}}({\bm{X}},A))(\mu({\bm{X}},A)-\tilde{\mu}({\bm{X}},A))]\rVert+O_{p}(n^{-1/2}),$

where (a) is by the property of the targeted regularization, namely,
$\frac{1}{n}\sum_{i=1}^{n}\hat{{\bm{w}}}({\bm{X}}_{i},A_{i})(Y_{i}-\tilde{\mu}({\bm{X}}_{i},A_{i}))=0$;
(b) is because of the uniform concentration of the empirical process, again
using the vanishing Rademacher complexity and uniform boundedness; (c)
integrates over $y$; (d) uses the definition of $\psi$ and the importance
sampling formula with ${\bm{w}}$. Since the link function $g$ is continuously
differentiable, invertible and strictly monotone, then by the mean value
theorem there exists ${\bm{\epsilon}}^{\prime}\in(0,\hat{{\bm{\epsilon}}})$
such that

$\displaystyle\tilde{\mu}({\bm{X}},A)$
$\displaystyle=g^{-1}(g(\hat{\mu}({\bm{X}},A))+\hat{{\bm{\epsilon}}})=\hat{\mu}({\bm{X}},A)+(g^{-1})^{\prime}(g(\hat{\mu}({\bm{X}},A))+{\bm{\epsilon}}^{\prime})\hat{{\bm{\epsilon}}}.$

From the uniform boundedness and smoothness of the link function, we have that
$\hat{c}=(g^{-1})^{\prime}(g(\hat{\mu}({\bm{X}},A)+{\bm{\epsilon}}^{\prime})<C$
for some constant $C>0$. Then, using the above result in the last term of
Equation 19, we obtain

$\displaystyle\lVert\mathbb{E}[(\hat{{\bm{w}}}({\bm{X}},A)-{\bm{w}}({\bm{X}},A))(\mu({\bm{X}},A)-\tilde{\mu}({\bm{X}},A))]\rVert$
(20)
$\displaystyle\quad\leq\mathbb{E}[(\hat{{\bm{w}}}({\bm{X}},A)-{\bm{w}}({\bm{X}},A))(\mu({\bm{X}},A)-\hat{\mu}({\bm{X}},A))]+C\lVert(\hat{{\bm{w}}}({\bm{X}},A)-{\bm{w}}({\bm{X}},A))\hat{{\bm{\epsilon}}}\rVert$
$\displaystyle\quad\leq
O_{p}(r_{1}(n)r_{2}(n))+O_{p}(r_{2}(n))\lVert\hat{{\bm{\epsilon}}}\rVert$

To complete the proof, we will show that
$\lVert\hat{{\bm{\epsilon}}}\rVert=O_{p}(r_{1}(n))+O_{p}(n^{-1/2})$. Letting
$\hat{c}_{i}$ be as in the Taylor expansion above, we can re-arrange the
targeted regularization condition such that

${\bm{0}}=\frac{d}{d{\bm{\epsilon}}}{\mathcal{R}}^{\text{tr}}(\hat{\mu},\hat{{\bm{w}}},\hat{{\bm{\epsilon}}})=\frac{1}{n}\sum_{i=1}^{n}\hat{{\bm{w}}}({\bm{X}},A)(\hat{\mu}({\bm{X}},A)-Y)+\frac{1}{n}\sum_{i=1}^{n}\hat{{\bm{w}}}({\bm{X}}_{i},A_{i})\hat{c}_{i}\hat{{\bm{\epsilon}}}.$

Hence, we can write $\hat{{\bm{\epsilon}}}$ with the closed-form expression

$\displaystyle\hat{{\bm{\epsilon}}}$
$\displaystyle=\operatorname*{arg\,min}_{{\bm{\epsilon}}}{\mathcal{R}}^{\text{tr}}(\hat{\mu},\hat{{\bm{w}}},{\bm{\epsilon}})=\frac{n^{-1}\sum_{i=1}^{n}\hat{{\bm{w}}}({\bm{X}}_{i},A_{i})(Y_{i}-\hat{\mu}({\bm{X}}_{i},A_{i}))}{n^{-1}\sum_{i=1}^{n}\hat{c}_{i}\hat{{\bm{w}}}({\bm{X}}_{i},A_{i})^{2}}.$

Since the denominator is uniformly bounded in a neighborhood of the solution
as in $\lVert 1/{\mathcal{W}}\rVert_{\infty}<\infty$, and $g$ is strictly
monotone and continuously differentiable, implies that $\hat{c}_{i}$ is
uniformly lower bounded. Hence, there is $C^{\prime}>0$ such that

$\displaystyle\lVert\hat{{\bm{\epsilon}}}\rVert$ $\displaystyle\leq
C^{\prime}\lVert\textstyle{n^{-1}\sum_{i=1}^{n}\hat{{\bm{w}}}({\bm{X}}_{i},A_{i})(Y_{i}-\hat{\mu}({\bm{X}}_{i},A_{i}))}\rVert$
(21)
$\displaystyle\leq^{(a)}C^{\prime}\lVert\mathbb{E}[\hat{{\bm{w}}}({\bm{X}},A)(Y-\hat{\mu}({\bm{X}},A))\rVert+O_{p}(n^{-1/2})$
$\displaystyle\leq^{(b)}O_{p}(r_{1}(n))+O_{p}(n^{-1/2}),$

where (a) uses the uniform concentration of the empirical process and (b)
again uses the uniform boundedness of $\hat{{\bm{w}}}$. The proof now follows
from combining equations (19), (20), and (21). ∎

##  Appendix C Another Example of SRF vs ERF

The following examples show a simple case where the SRF and ERF estimands are
different, thereby demonstrating why the ERF is not a useful estimate of the
effect of an exposure shift. Consider a setting in which $\mu({\bm{X}},A)=AX$
with $X\sim N(0,1)$, $A\sim N(X,1)$. Now consider the exposure shift induced
by $\tilde{A}=cA$ for some $c\in\mathbb{R}$. Using the SRF formulation, we
find that
${\bm{\psi}}=\mathbb{E}[\mu({\bm{X}},\tilde{A})]=\mathbb{E}[\mathbb{E}[{\bm{X}}(cA)|{\bm{X}}]]=c\mathbb{E}[{\bm{X}}^{2}]=c$.
On the other hand, the ERF is
$\xi(a)=\mathbb{E}[\mu({\bm{X}},A)|A=a]=a\mathbb{E}[{\bm{X}}]=0$ for all
$a\in\mathbb{R}$. Therefore, estimators of the two estimands return two
different estimates. Moreover, the ERF is identically zero for every treatment
value. Thus, it cannot be used to approximate the value of the effect of the
exposure shift, even when ${\bm{\psi}}$ is correctly specified.

##  Appendix D Hardware/Software/Data Access

We ran all of our experiments both in the simulation study and the application
section using Pytorch (Paszke et al., 2019) on a high-performance computing
cluster equipped with Intel 8268 ”Cascade Lake” processors. Due to the
relatively small size of the datasets, hardware limitations, and the large
number of simulations required, we did not require the use of GPUs. Instead we
found that using only CPUs run in parallel sufficed. Reproducing the full set
of experiments takes approximately 12 hours with 100 parallel processes, each
with 4 CPU cores and 8GB of RAM. Each process runs a different random seed for
the experiment configuration.

The code for reproducibility is provided on the submission repository along
with the data sources for the simulation experiments. The datasets for these
experiments were obtained from the public domain and were adapted from the
GitHub repositories shared by Nie et al. (2021) and Bica et al. (2020) as
explained in the experiment details section. The data for the application was
purchased from https://resdac.org/. Due to a data usage agreement and privacy
concerns, manipulation of these data requires IRB approval under which the
authors have completed the training and for which reason the data cannot be
shared with the public.