On the Evaluation of Surrogate Markers in Real World Data Settings–Double
Robustness April

10.1111/j.1541-0420.2005.00454.x

# On the Evaluation of Surrogate Markers in Real World Data Settings

Larry Han∗<EMAIL_ADDRESS>Xuan Wang∗∗<EMAIL_ADDRESS>and
Tianxi Cai∗∗∗
Department of Biostatistics<EMAIL_ADDRESS>Harvard T.H. Chan School of
Public Health  Boston  Massachusetts  02115

(Received April 2021\. Revised August 2021\. Accepted December 2021.; 2021)

###### Abstract

Shortcomings of randomized clinical trials are pronounced in urgent health
crises, when rapid identification of effective treatments is critical.
Leveraging short-term surrogates in real-world data (RWD) can guide
policymakers evaluating new treatments. In this paper, we develop novel
estimators for the proportion of treatment effect (PTE) on the true outcome
explained by a surrogate in RWD settings. We propose inverse probability
weighted and doubly robust (DR) estimators of an optimal transformation of the
surrogate and PTE by semi-nonparametrically modeling the relationship between
the true outcome and surrogate given baseline covariates. We show that our
estimators are consistent and asymptotically normal, and the DR estimator is
consistent when either the propensity score model or outcome regression model
is correctly specified. We compare our proposed estimators to existing
estimators and show a reduction in bias and gains in efficiency through
simulations. We illustrate the utility of our method in obtaining an
interpretable PTE by conducting a cross-trial comparison of two biologic
therapies for ulcerative colitis.

###### keywords:

Double robustness; Proportion of treatment effect explained; Real world data;
Semi-nonparametric estimation; Surrogate marker.

††volume: 64

##  1 Introduction

While randomized clinical trials (RCTs) remain the gold standard instrument
for identifying efficacious and safe drugs (Concato, et al., 2000), RCTs
typically require long-term follow-up of patients to observe a sufficient
number of events to estimate treatment effects and even then, the outcomes may
be costly to measure (Bentley et al., 2019). RCTs are also often limited to
narrowly defined patient populations with results that are not always
generalizable. These shortcomings are especially pronounced in urgent health
crises and have led to increased interest in using real world data (RWD) and
shorter term surrogate endpoints to efficiently and effectively evaluate
treatments (Hernán & Robins, 2016; Corrigan-Curay et al., 2018; Hey et al.,
2020; Gyawali et al., 2020). For example, as multiple pharmaceutical companies
pushed to develop vaccines and treatments for COVID-19 and to test them in
RCTs (Lurie et al., 2020), RWD from electronic health records (EHRs) have been
collected at breakneck speed (Brat et al., 2020).

The use of valid surrogate markers to infer treatment effects on long term
outcomes has the potential to reduce trial cost and study duration (Ciani et
al., 2017; Wickström & Moseley, 2017). The explosion in recent years of RWD
highlights an untapped opportunity to identify and validate surrogate markers.
Since Prentice (1989) originally proposed a definition and operational
criteria for identifying valid surrogate markers, many statistical methods
have been developed to make inference about the proportion of treatment effect
(PTE) explained by a surrogate in RCT settings (Freedman & Schatzkin, 1992;
Lin et al., 1997; Wang & Taylor, 2002; Parast et al., 2016; Price et al.,
2018; Wang et al., 2020). For example, Freedman & Schatzkin (1992) proposed a
parametric model-based estimate assuming two regression models which rarely
hold simultaneously (Lin et al., 1997). Wang & Taylor (2002) proposed
alternative measures of PTE that examined what the treatment effect would have
been if the surrogate had the same distribution across treatment groups.
Parast et al. (2016) proposed a fully nonparametric estimation procedure for
the PTE defined in Wang & Taylor (2002). More recently, Wang et al. (2020)
proposed an alternative non-parametric PTE estimator by identifying an optimal
transformation of the surrogate $S$, $g_{\scriptscriptstyle\sf opt}:S\to
g_{\scriptscriptstyle\sf opt}(S)$, such that $g_{\scriptscriptstyle\sf
opt}(S)$ optimally predicts $Y$. In addition to requiring weaker assumptions
than those required by Parast et al. (2016), this approach has the advantage
of providing a direct approximation to the treatment effect on the outcome $Y$
using the treatment effect on $g_{\scriptscriptstyle\sf opt}(S)$.

These existing PTE estimates are derived for data from RCTs and not directly
applicable to RWD where treatment assignment $A$ may depend on confounding
factors $\mathbf{X}$. In this paper, we follow the strategy of Wang et al.
(2020) to define PTE based on $g_{\scriptscriptstyle\sf opt}$ and propose both
IPW and doubly robust PTE estimators using RWD by semi-non-parametrically
modeling the relationship between $Y$ and $S$ given $\mathbf{X}$ and imposing
a propensity score (PS) model for $A\mbox{ given }\mathbf{X}$. We propose
perturbation resampling methods for variance and confidence interval
estimation. We establish the asymptotic properties of the proposed estimators,
including double robustness of the proposed estimator in that it is consistent
when either the PS model or the outcome regression (OR) models is correctly
specified. Our simulation studies demonstrate that the proposed estimators and
inference procedures perform well in finite samples. Finally, we illustrate
the utility of our proposed procedures by conducting a cross-trial comparison
of two biologic therapies for ulcerative colitis (UC).

##  2 Methods

###  2.1 Setting and notations

Let $Y$ be the primary outcome and $S$ be the surrogate marker, both of which
may be discrete or continuous. Throughout, the notation takes $S$ to be
continuous, but all derivations and theoretical results remain valid if $S$ is
discrete by replacing density functions with probability mass functions. We
denote $\\{Y^{(a)},S^{(a)}\\}$ as the respective potential primary outcome and
surrogate marker under treatment $A=a$, where $A=1$ and $A=0$ denote the
treatment and the control group, respectively. With RWD, only
$Y_{i}=A_{i}Y_{i}^{(1)}+(1-A_{i})Y_{i}^{(0)}$ and
$S_{i}=A_{i}S_{i}^{(1)}+(1-A_{i})S_{i}^{(0)}$ can be observed for an
individual $i$, and the treatment assignment $A_{i}$ may depend on baseline
confounding factors $\mathbf{X}_{i}$. For identifiability, we require the
standard assumptions (Rubin, 2005; Imbens & Rubin, 2015):

${\pi_{a}(\mathbf{x})\equiv P(A=a|\mathbf{X}=\mathbf{x})\in(0,1)}$ (1)
$\left(Y^{(1)},Y^{(0)},S^{(1)},S^{(0)}\right)\perp A\mid\mathbf{X}$ (2)

Assumption (1) states that within all covariate levels, patients may receive
either treatment so that the PS is bounded away from $0$ and $1$. Assumption
(2) implies that $\mathbf{X}$ includes all confounders that can affect the
primary outcome and treatment simultaneously, or the surrogate and treatment
simultaneously (Rubin, 2005; Imbens & Rubin, 2015). We assume that the RWD for
analysis consist of $n$ independent and identically distributed random
variables
$\\{\mathbf{D_{i}}=(Y_{i},S_{i},A_{i},\mathbf{X}_{i}^{\scriptscriptstyle\sf
T})^{\scriptscriptstyle\sf T},i=1,...,n\\}$.

###  2.2 Target parameter and leveraging surrogates

The average treatment effect on $Y$ is defined as:

$\Delta=\mu_{1}-\mu_{0},\quad\text{ where }\mu_{a}=E(Y^{(a)})=\int E(Y\mid
A=a,\mathbf{X})d\mathbb{F}(\mathbf{X}),$

and $\mathbb{F}(\mathbf{x})=P(\mathbf{X}\leq\mathbf{x})$. Without loss of
generality, we assume that $\Delta>0$. To approximate $\Delta$ based on the
treatment effect on $S$, Wang et al. (2020) identified a transformation
function $g_{\scriptscriptstyle\sf opt}(\cdot)$ such that the treatment effect
on the transformed surrogate, $\Delta_{g_{\scriptscriptstyle\sf
opt}}=E[g_{opt}(S^{(1)})-g_{opt}(S^{(0)})],$ can optimally predict $\Delta$ in
a certain sense. More formally, the optimality of $g_{\scriptscriptstyle\sf
opt}$ is with respect to minimizing the mean squared error

$\displaystyle\mathcal{L}_{oracle}(g_{opt})=E\left[\left(Y^{(1)}-Y^{(0)}\right)-\left\\{g_{{\text{opt}}}\left(S^{(1)}\right)-g_{\text{opt}}\left(S^{(0)}\right)\right\\}\right]^{2}$

under the working assumption of $(Y^{(1)},S^{(1)})\perp(Y^{(0)},S^{(0)})$. It
was shown that $g_{opt}$ takes the form

$g_{opt}(s)=m(s)+\lambda\mathcal{P}_{0}(s)\quad\mbox{with}\quad
m(s)=m_{1}(s)\mathcal{P}_{1}(s)+m_{0}(s)\mathcal{P}_{0}(s),$

where $m_{a}(s)=E(Y^{(a)}|S^{(a)}=s)$,
$\mathcal{P}_{a}(s)={f_{a}(s)}(f_{0}(s)+f_{1}(s))^{-1}$,
$f_{a}(s)={dF_{a}(s)}/{ds},$

$\
\lambda=\frac{\int\\{m_{0}(s)-m_{1}(s)\\}\mathcal{P}_{1}(s)dF_{0}(s)}{\int\mathcal{P}_{0}(s)dF_{0}(s)}=\frac{\mu_{0}-\int
m(s)dF_{0}(s)}{\int\mathcal{P}_{0}(s)dF_{0}(s)},$

and $F_{a}(s)=P(S^{(a)}\leq s)$. In addition, by employing the transformation
$g_{\scriptscriptstyle\sf opt}$ and defining the PTE of $S$ as
$\mbox{PTE}_{g_{\scriptscriptstyle\sf
opt}}\equiv{\Delta_{g_{\scriptscriptstyle\sf opt}}}/{\Delta},$ Wang et al.
(2020) showed that even if the working independence assumption does not hold,
$\mbox{PTE}_{g_{\scriptscriptstyle\sf opt}}\in[0,1]$ provided that

$\displaystyle\text{(A1)}\quad\mathbb{S}_{1}(u)\geq\mathbb{S}_{0}(u)\quad\mbox{for
all $u$},$
$\displaystyle\text{(A2)}\quad\mathbb{M}_{1}(u)\geq\mathbb{M}_{0}(u)\quad\mbox{for
all $u$ in the common support of $g_{opt}(S^{(1)})$ and $g_{opt}(S^{(0)})$},$

where $\mathbb{S}_{a}(u)=P\\{g_{opt}(S^{(a)})\geq u\\}$ and
$\mathbb{M}_{a}(u)=E(Y^{(a)}|g_{opt}(S^{(a)})=u)$, for $a=0,1$. Assumptions
(A1) and (A2) are weaker than those required in the literature to ensure that
the PTE is between 0 and 1 and hence to avoid the surrogate paradox
(VanderWeele, 2013). Our goal here is to construct robust estimates for
$g_{\scriptscriptstyle\sf opt}$ and $\mbox{PTE}_{g_{\scriptscriptstyle\sf
opt}}$ using RWD.

##  3 Two Proposed Estimation Methods

Estimation of $g_{opt}(s)$ and $\mbox{PTE}_{g_{\scriptscriptstyle\sf opt}}$
using RWD is more challenging than using RCT data because we cannot directly
estimate $m(s)$, $\lambda$, and $\mathcal{P}_{a}(s)$ due to confounding. We
propose an inverse probability weighted (IPW) estimator and a doubly robust
(DR) estimator for $g_{\scriptscriptstyle\sf opt}$ and
$\mbox{PTE}_{g_{\scriptscriptstyle\sf opt}}$ accounting for the effects of
$\mathbf{X}$ on $A$, $Y$ and $S$. We first present the simpler IPW estimator
and then the DR estimator. For both estimators, we impose a parametric model
for $\pi_{1}(\mathbf{X})$, denoted by
$\pi_{1}(\mathbf{X};\boldsymbol{\alpha})$, where $\boldsymbol{\alpha}$ is a
finite dimensional parameter that can be estimated as the standard maximum
likelihood estimator, $\widehat{\boldsymbol{\alpha}}$. A simple example is a
logistic regression model
$\pi_{1}(\mathbf{X};\boldsymbol{\alpha})=G\\{\boldsymbol{\alpha}^{\scriptscriptstyle\sf
T}\boldsymbol{\Phi}(\mathbf{X})\\}$, where $G(x)=e^{x}/(1+e^{x})$ and
$\boldsymbol{\Phi}(\mathbf{X})$ is a vector of basis functions of $\mathbf{X}$
to account for potential non-linear effects.

###  3.1 IPW Estimation

To construct an IPW estimator for $g_{\scriptscriptstyle\sf opt}$, we first
obtain IPW kernel smoothed estimators for $m_{a}(s)$ and $f_{a}(s)$
respectively as

$\widehat{m}_{a}(s)=\frac{\sum_{i=1}^{n}K_{h}(S_{i}-s)Y_{i}\widehat{\omega}_{ai}}{\sum_{i=1}^{n}K_{h}(S_{i}-s)\widehat{\omega}_{ai}}\quad\mbox{and}\quad\widehat{f}_{a}(s)=\frac{\sum_{i=1}^{n}K_{h}(S_{i}-s)\widehat{\omega}_{ai}}{\sum_{i=1}^{n}\widehat{\omega}_{ai}},$

where
$\widehat{\omega}_{ai}=I(A_{i}=a)/\pi_{a}(\mathbf{X}_{i},\widehat{\boldsymbol{\alpha}}),$
$K_{h}(\cdot)=h^{-1}K(\cdot/h)$, $K(\cdot)$ is a symmetric density function
and $h=O(n^{-\nu})$ with $\nu\in(1/4,1/2)$. Then $m(\cdot)$,
$\mathcal{P}_{a}(\cdot)$ and $\lambda$ may be estimated as

$\widehat{m}(s)=\sum_{a=0}^{1}\widehat{m}_{a}(s)\widehat{\mathcal{P}}_{a}(s),\quad\widehat{\mathcal{P}}_{a}(s)=\frac{\widehat{f}_{a}(s)}{\widehat{f}_{1}(s)+\widehat{f}_{0}(s)},\quad\hat{\lambda}=\frac{\int(\widehat{m}_{0}(s)-\widehat{m}_{1}(s))\widehat{\mathcal{P}}_{1}(s)\widehat{f}_{0}(s)ds}{\int\widehat{\mathcal{P}}_{0}(s)\widehat{f}_{0}(s)ds},$

respectively. Subsequently, we construct plug-in estimators for
$g_{\scriptscriptstyle\sf opt}(s)$, $\Delta_{g_{\scriptscriptstyle\sf opt}}$,
$\Delta$ and $\mbox{PTE}_{g_{\scriptscriptstyle\sf opt}}$ as

$\widehat{g}(s)=\widehat{m}(s)+\hat{\lambda}\widehat{\mathcal{P}}_{0}(s),\quad\widehat{\Delta}_{\widehat{g}}=\widehat{\mu}_{1,\widehat{g}}-\widehat{\mu}_{0,\widehat{g}},\quad\widehat{\Delta}=\widehat{\mu}_{1}-\widehat{\mu}_{0},\quad\mbox{and}\quad\widehat{\mbox{PTE}}_{\widehat{g}}=\frac{\widehat{\Delta}_{\widehat{g}}}{\widehat{\Delta}},$

where

$\widehat{\mu}_{a,g}=\frac{\sum_{i=1}^{n}g(S_{i})\widehat{\omega}_{ai}}{\sum_{i=1}^{n}\widehat{\omega}_{ai}}\quad\mbox{and}\quad\widehat{\mu}_{a}=\frac{\sum_{i=1}^{n}Y_{i}\widehat{\omega}_{ai}}{\sum_{i=1}^{n}\widehat{\omega}_{ai}.}$

We show in Appendix 1 of the supplementary materials that when
$\pi_{1}(\mathbf{x};\boldsymbol{\alpha})$ is correctly specified,
$\widehat{\mbox{PTE}}_{\widehat{g}}$ is consistent for
$\mbox{PTE}_{g_{\scriptscriptstyle\sf opt}}$. We also show that
$\sqrt{n}(\widehat{\mbox{PTE}}_{\widehat{g}}-\mbox{PTE}_{g_{\scriptscriptstyle\sf
opt}})$ converges in distribution to a normal distribution with mean $0$ and
variance $\sigma^{2}$, where the form of $\sigma^{2}$ is given in the
supplementary materials.

###  3.2 Doubly Robust Estimation

When the PS model is incorrectly specified, the IPW estimator is likely to be
biased. Here, we propose augmented IPW estimators for
$g_{\scriptscriptstyle\sf opt}$ and $\mbox{PTE}_{g_{\scriptscriptstyle\sf
opt}}$ to achieve improved robustness and efficiency. Following Robins et al.
(1994), for any counterfactual random variable $U^{(a)}$, an augmented IPW
estimator for its mean $E(U^{(a)})$ can be constructed as

$n^{-1}\sum_{i=1}^{n}\left\\{\widehat{\omega}_{ai}U_{i}-(\widehat{\omega}_{ai}-1)\widehat{\phi}_{a}(\mathbf{X}_{i})\right\\},$

where $\widehat{\phi}_{a}(\mathbf{X}_{i})$ is an estimator for
$E(U_{i}^{(a)}\mid\mathbf{X}_{i})$ derived under a specified model. This
estimator is doubly robust in the sense that it is consistent for $E(U^{(a)})$
when either the PS model for $\pi_{a}(\mathbf{X})$ or the outcome model for
$E(U_{i}^{(a)}\mid\mathbf{X}_{i})$ is correctly specified. Deriving an
augmented IPW estimator for $\mbox{PTE}_{g_{\scriptscriptstyle\sf opt}}$ is
more involved since $g_{\scriptscriptstyle\sf opt}(S)$ involves conditional
mean functions of $Y^{(a)}\mid S^{(a)}$ and density functions of $S^{(a)}$ for
$a=0,1$.

To construct a DR estimator for $g_{\scriptscriptstyle\sf
opt}(s)=m(s)+\lambda\mathcal{P}_{0}(s)$, we propose the following DR
estimators for $m_{a}(s)$ and $f_{a}(s)$ respectively,

$\displaystyle\widehat{m}_{\scriptscriptstyle a,\sf DR}(s)$
$\displaystyle=\frac{\widehat{\mathcal{M}}_{\scriptscriptstyle a,\sf
DR}(s)}{\widehat{f}_{\scriptscriptstyle a,\sf DR}(s)},$ (3a)
$\displaystyle\widehat{\mathcal{M}}_{\scriptscriptstyle a,\sf DR}(s)$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\left\\{K_{h}(S_{i}-s)Y_{i}\widehat{\omega}_{ai}-(\widehat{\omega}_{ai}-1)\widehat{\psi}^{\scriptscriptstyle{\dagger}}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{X}_{i})\widehat{\psi}^{\scriptscriptstyle{\dagger}}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{X}_{i})\right\\},$ (3b)
$\displaystyle\widehat{f}_{\scriptscriptstyle a,\sf DR}(s)$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\left\\{K_{h}(S_{i}-s)\widehat{\omega}_{ai}-(\widehat{\omega}_{ai}-1)\widehat{\psi}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{X}_{i})\right\\},$ (3c)

where $\widehat{\psi}_{a,\scriptscriptstyle\sf m}(\mathbf{x})$ and
$\widehat{\psi}_{a,\scriptscriptstyle\sf f}(s;\mathbf{x})$ are the respective
estimators for

$\displaystyle\psi_{a,\scriptscriptstyle\sf m}(s;\mathbf{x})$
$\displaystyle=E(Y_{i}^{(a)}\mid
S_{i}^{(a)}=s,\mathbf{X}_{i}=\mathbf{x})=E(Y_{i}\mid
A_{i}=a,S_{i}=s,\mathbf{X}_{i}=\mathbf{x})\ \mbox{and }$
$\displaystyle\psi_{a,\scriptscriptstyle\sf f}(s;\mathbf{x})$
$\displaystyle=\frac{\partial P(S_{i}^{(a)}\leq
s\mid\mathbf{X}_{i}=\mathbf{x})}{\partial s},$

In Appendix 2 of the supplementary materials, we show that
$\widehat{m}_{\scriptscriptstyle a,\sf DR}(s)$ and
$\widehat{f}_{\scriptscriptstyle a,\sf DR}(s)$ are consistent for $m_{a}(s)$
and $f_{a}(s)$ if either
$\sup_{\mathbf{x}}|\widehat{\pi}_{a}(\mathbf{x})-\pi_{a}(\mathbf{x})|\to 0$ in
probability or $\sup_{\mathbf{x},s}\\{|\widehat{\psi}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{x})-\psi_{a,\scriptscriptstyle\sf
m}(s;\mathbf{x})|+|\widehat{\psi}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{x})-\psi_{a,\scriptscriptstyle\sf f}(s;\mathbf{x})|\\}\to 0$ in
probability.

To construct estimators $\widehat{\psi}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{x})$ and $\widehat{\psi}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{x})$, we impose flexible semi-non-parametric models for
$Y^{(a)}\mid S^{(a)},\mathbf{X}$ and $S^{(a)}\mid\mathbf{X}$ to minimize
assumptions on the dependency structure between $S$ and $Y$. Specifically, we
first impose a generalized regression model (GRM) (Han, 1987) for $S_{i}\mid
A_{i}=a,\mathbf{X}_{i}$:

$S_{i}=\mathcal{D}_{a}\odot\mathcal{H}_{a}(\mathbf{X}_{i}^{\scriptscriptstyle\sf
T}\boldsymbol{\gamma}_{a},\epsilon_{ia})\quad\mbox{with}\quad
P(\epsilon_{ai}\leq e\mid\mathbf{X}_{i})=\mathcal{F}_{a}(e),$

where $\mathcal{D}_{a}(\cdot)$ is an increasing function and
$\mathcal{H}_{a}(\cdot,\cdot)$ is a strictly increasing function of each of
its arguments, and the unknown covariate effects
$\boldsymbol{\gamma}_{a}=(\gamma_{a1},...,\gamma_{ap})^{\scriptscriptstyle\sf
T}$ are constrained to the unit sphere
$\Omega:\\{\boldsymbol{\gamma}:\|\boldsymbol{\gamma}\|_{2}=1\\}$ for
identifiability. With the given $\boldsymbol{\gamma}_{a}$ under GRM and the
no-unmeasured-confounders assumption, $\psi_{a,\scriptscriptstyle\sf
f}(s;\mathbf{x})$ can be estimated non-parametrically via kernel smoothing. To
estimate $\boldsymbol{\gamma}_{a}$, Sherman (1993) showed that the maximum
rank correlation estimator

$\widehat{\boldsymbol{\gamma}}_{a}=\mbox{argmax}_{\boldsymbol{\gamma}\in\Omega}\left\\{\sum_{i\neq
j,A_{i}=A_{j}=a}I(\mathbf{X}_{i}^{\scriptscriptstyle\sf
T}\boldsymbol{\gamma}>\mathbf{X}_{j}^{\scriptscriptstyle\sf
T}\boldsymbol{\gamma})I(S_{i}>S_{j})\right\\}$

is consistent and asymptotically normal for $\boldsymbol{\gamma}_{a}$.
Subsequently, we estimate $\psi_{a,\scriptscriptstyle\sf f}(s,\mathbf{x})$ as

$\widehat{\psi}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{x})=\frac{\sum_{i=1}^{n}K_{\zeta}(\widehat{\boldsymbol{\gamma}}_{a}\mathbf{X}_{i}-\widehat{\boldsymbol{\gamma}}_{a}^{\scriptscriptstyle\sf
T}\mathbf{x})K_{h}(S_{i}-s)}{\sum_{i=1}^{n}K_{\zeta}(\widehat{\boldsymbol{\gamma}}_{a}\mathbf{X}_{i}-\widehat{\boldsymbol{\gamma}}_{a}^{\scriptscriptstyle\sf
T}\mathbf{x})}.$ (4)

To estimate $\psi_{a,\scriptscriptstyle\sf m}(s;\mathbf{x})$, we impose a
varying-coefficient generalized linear model (VGLM) (Hastie & Tibshirani,
1993):

$E(Y_{i}\ \mid
A_{i}=a,S_{i}=s,\mathbf{X}_{i})=M\\{\boldsymbol{\beta}_{a}(S_{i})^{\scriptscriptstyle\sf
T}\vec{\mathbf{X}}_{i}\\},$

where $M(\cdot)$ is a known smooth link function,
$\vec{\mathbf{x}}=(1,\mathbf{x}^{\scriptscriptstyle\sf
T})^{\scriptscriptstyle\sf T}$ for any vector $\mathbf{x}$ and
$\boldsymbol{\beta}_{a}(s)$ is an unknown $p+1$ dimensional unspecified smooth
coefficient functions. We may estimate $\boldsymbol{\beta}_{a}(s)$ as
$\widehat{\boldsymbol{\beta}}_{a}(s)$, the solution to

$\widehat{\mathbf{U}}_{a}(\boldsymbol{\beta};s)=n^{-1}\sum_{i=1}^{n}I(A_{i}=a)K_{h}(S_{i}-s)\vec{\mathbf{X}}_{i}\left\\{Y_{i}-M(\boldsymbol{\beta}^{\scriptscriptstyle\sf
T}\vec{\mathbf{X}}_{i})\right\\}=\mathbf{0}$

Then we estimate $\psi_{a,\scriptscriptstyle\sf m}(s;\mathbf{x})$ as

$\widehat{\psi}_{a,\scriptscriptstyle\sf
m}(s,\mathbf{x})=M\\{\widehat{\boldsymbol{\beta}}_{a}(s)^{\scriptscriptstyle\sf
T}\vec{\mathbf{x}}\\}.$ (5)

These estimators (4) and (5) can then be plugged into (3b) and (3c) to
construct $\widehat{f}_{\scriptscriptstyle a,\sf DR}(s)$ and
$\widehat{m}_{\scriptscriptstyle a,\sf DR}(s)$ as in (3a).

Based on $\widehat{m}_{\scriptscriptstyle a,\sf DR}(s)$ and
$\widehat{f}_{\scriptscriptstyle a,\sf DR}(s)$, we obtain a doubly robust
estimator for $g_{opt}(s)$ as:

$\widehat{g}_{\scriptscriptstyle\sf DR}(s)=\widehat{m}_{\scriptscriptstyle\sf
DR}(s)+\widehat{\lambda}_{\scriptscriptstyle\sf
DR}\widehat{\mathcal{P}}_{\scriptscriptstyle 0,\sf DR}(s),$

where $\widehat{m}_{\scriptscriptstyle\sf
DR}(s)=\sum_{a=0}^{1}\widehat{m}_{\scriptscriptstyle a,\sf
DR}(s)\widehat{\mathcal{P}}_{\scriptscriptstyle a,\sf DR}(s)$,

$\widehat{\lambda}_{\scriptscriptstyle\sf
DR}=\frac{\int\left\\{\widehat{m}_{\scriptscriptstyle 0,\sf
DR}(s)-\widehat{m}_{\scriptscriptstyle 1,\sf
DR}(s)\right\\}\widehat{\mathcal{P}}_{\scriptscriptstyle 1,\sf
DR}(s)\widehat{f}_{\scriptscriptstyle 0,\sf
DR}(s)ds}{\int\widehat{\mathcal{P}}_{\scriptscriptstyle 0,\sf
DR}(s)\widehat{f}_{\scriptscriptstyle 0,\sf
DR}(s)ds},\quad\mbox{and}\quad\widehat{\mathcal{P}}_{\scriptscriptstyle a,\sf
DR}(s)=\frac{\widehat{f}_{\scriptscriptstyle a,\sf
DR}(s)}{\widehat{f}_{\scriptscriptstyle 0,\sf
DR}(s)+\widehat{f}_{\scriptscriptstyle 1,\sf DR}(s)},$

for $a=0,1$. We can now construct a plug-in estimator for
$\Delta_{g_{\scriptscriptstyle\sf opt}}$ as
$\widehat{\Delta}_{\widehat{g},\scriptscriptstyle\sf
DR}=\widehat{\mu}_{1,\widehat{g},\scriptscriptstyle\sf
DR}-\widehat{\mu}_{0,\widehat{g},\scriptscriptstyle\sf DR},$ where

$\widehat{\mu}_{a,g,\scriptscriptstyle\sf
DR}=n^{-1}\sum_{i=1}^{n}\left\\{g(S_{i})\widehat{\omega}_{ai}-(\widehat{\omega}_{ai}-1)\hat{\zeta}_{a,g}(\mathbf{X}_{i})\right\\},$

and $\hat{\zeta}_{a,g}(\mathbf{x})=\int
g(s)\widehat{\psi}_{a,\scriptscriptstyle\sf f}(s,\mathbf{x})ds$ is an
estimator for
$\zeta_{a,g}(\mathbf{x})=E\\{g(S_{i}^{(a)})\mid\mathbf{X}_{i}=\mathbf{x}\\}$
derived under the GRM. Similarly, we obtain
$\widehat{\Delta}_{\scriptscriptstyle\sf DR}=\widehat{\mu}_{\scriptscriptstyle
1,\sf DR}-\widehat{\mu}_{\scriptscriptstyle 0,\sf DR}$ to estimate $\Delta$,
where

$\widehat{\mu}_{\scriptscriptstyle a,\sf
DR}=n^{-1}\sum_{i=1}^{n}\left\\{Y_{i}\widehat{\omega}_{ai}-(\widehat{\omega}_{ai}-1)\hat{\zeta}_{a}(\mathbf{X}_{i})\right\\},$

where $\hat{\zeta}_{a}(\mathbf{x})=\int\widehat{\psi}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{x})\widehat{\psi}_{a,\scriptscriptstyle\sf f}(s;\mathbf{x})ds$ is
an estimator for
$\zeta_{a}(\mathbf{x})=E(Y_{i}^{(a)}\mid\mathbf{X}_{i}=\mathbf{x})$. Finally,
we estimate $\mbox{PTE}_{g_{\scriptscriptstyle\sf opt}}$ as
$\widehat{\mbox{PTE}}_{\widehat{g},\scriptscriptstyle\sf
DR}={\widehat{\Delta}_{\widehat{g},\scriptscriptstyle\sf
DR}}/{\widehat{\Delta}}_{\scriptscriptstyle\sf DR}.$ Following similar
arguments as given in Appendix 2 of the supplementary materials, it is not
difficult to show that $\widehat{\Delta}_{\widehat{g},\scriptscriptstyle\sf
DR}$, $\widehat{\Delta}_{\scriptscriptstyle\sf DR}$ and
$\widehat{\mbox{PTE}}_{\widehat{g},\scriptscriptstyle\sf DR}$ are doubly
robust estimators for $\Delta_{g_{\scriptscriptstyle\sf opt}}$, $\Delta$, and
$\mbox{PTE}_{g_{\scriptscriptstyle\sf opt}}$, respectively.

##  4 Perturbation Resampling

We propose to estimate the variability and construct confidence intervals of
our proposed estimators using a perturbation-resampling approach (Jin et al.,
2001; Tian et al., 2005). For resampling, we generate
$\\{\mathbf{V}^{\scriptscriptstyle[b]}=(V_{1}^{\scriptscriptstyle[b]},...,V_{n}^{\scriptscriptstyle[b]})^{\scriptscriptstyle\sf
T},b=1,...,B\\}$, which are $n\times B$ independent and identically
distributed non-negative random variables from a known distribution with unit
mean and unit variance, such as the unit exponential distribution. For the IPW
estimators, for each set of
$\mathbf{V}=(V_{1},...,V_{n})^{\scriptscriptstyle\sf T}$, we let
$\bar{\mathcal{V}}_{i}=V_{i}/(n^{-1}\sum_{i=1}^{n}V_{i})$,

$\widehat{m}_{a}^{*}(s)=\frac{\sum_{i=1}^{n}K_{h}(S_{i}-s)Y_{i}\bar{\mathcal{V}}_{i}\widehat{\omega}_{ai}^{*}}{\sum_{i=1}^{n}K_{h}(S_{i}-s)\bar{\mathcal{V}}_{i}\widehat{\omega}_{ai}^{*}},\quad\widehat{f}_{a}^{*}(s)=\frac{\sum_{i=1}^{n}K_{h}(S_{i}-s)\bar{\mathcal{V}}_{i}\widehat{\omega}_{ai}^{*}}{\sum_{i=1}^{n}\bar{\mathcal{V}}_{i}\widehat{\omega}_{ai}^{*}},\quad\widehat{\omega}_{ai}^{*}=\frac{I(A_{i}=a)}{\pi(\mathbf{X}_{i},\widehat{\boldsymbol{\alpha}}^{*})},$

where $\widehat{\boldsymbol{\alpha}}^{*}$ is obtained by fitting a weighted
logistic regression $A_{i}\sim G\\{\boldsymbol{\alpha}^{\scriptscriptstyle\sf
T}\boldsymbol{\Phi}(\mathbf{X}_{i})\\}$ with weights
$\\{\bar{\mathcal{V}}_{i}\\}$. The perturbed counterparts of
$\widehat{m}(\cdot)$, $\hat{\mathcal{P}}_{a}(\cdot)$ and $\widehat{\lambda}$
are obtained as

$\widehat{m}^{*}(s)=\sum_{a=0}^{1}\widehat{m}_{a}^{*}(s)\widehat{\mathcal{P}}_{a}^{*}(s),\quad\widehat{\mathcal{P}}_{a}^{*}(s)=\frac{\widehat{f}_{a}^{*}(s)}{\widehat{f}_{1}^{*}(s)+\widehat{f}_{0}^{*}(s)},\quad\widehat{\lambda}^{*}=\frac{\int\\{\widehat{m}_{0}^{*}(s)-\widehat{m}_{1}^{*}(s)\\}\widehat{\mathcal{P}}_{1}^{*}(s)\widehat{f}_{0}^{*}(s)ds}{\int\widehat{\mathcal{P}}_{0}^{*}(s)\widehat{f}_{0}^{*}(s)ds},$

respectively. Subsequently, we construct the perturbed counterparts of
$\widehat{g}(s)$, $\widehat{\Delta}_{g_{\scriptscriptstyle\sf opt}}$,
$\widehat{\Delta}$ and $\widehat{\mbox{PTE}}$ as

$\widehat{g}^{*}(s)=\widehat{m}^{*}(s)+\hat{\lambda}^{*}\widehat{\mathcal{P}}_{0}^{*}(s)\quad\widehat{\Delta}^{*}_{\widehat{g}^{*}}=\widehat{\mu}^{*}_{1,\widehat{g}^{*}}-\widehat{\mu}^{*}_{0,\widehat{g}^{*}},\quad\widehat{\Delta}^{*}=\widehat{\mu}^{*}_{1}-\widehat{\mu}^{*}_{0},\quad\mbox{and}\quad\widehat{\mbox{PTE}}^{*}_{\widehat{g}^{*}}=\frac{\widehat{\Delta}^{*}_{\widehat{g}^{*}}}{\widehat{\Delta}^{*}},$

where

$\widehat{\mu}^{*}_{a,g}=\frac{\sum_{i=1}^{n}g(S_{i})\bar{\mathcal{V}}_{i}\widehat{\omega}^{*}_{ai}}{\sum_{i=1}^{n}\bar{\mathcal{V}}_{i}\widehat{\omega}^{*}_{ai}}\quad\mbox{and}\quad\widehat{\mu}^{*}_{a}=\frac{\sum_{i=1}^{n}Y_{i}\bar{\mathcal{V}}_{i}\widehat{\omega}_{ai}^{*}}{\sum_{i=1}^{n}\bar{\mathcal{V}}_{i}\widehat{\omega}_{ai}^{*}}.$

For the DR estimators, for each set of $\mathbf{V}$, we let

$\displaystyle\widehat{m}^{*}_{\scriptscriptstyle a,\sf DR}(s)$
$\displaystyle=\frac{\widehat{\mathcal{M}}^{*}_{\scriptscriptstyle a,\sf
DR}(s)}{\widehat{f}^{*}_{\scriptscriptstyle a,\sf DR}(s)},$
$\displaystyle\widehat{\mathcal{M}}^{*}_{\scriptscriptstyle a,\sf DR}(s)$
$\displaystyle=n^{-1}\sum_{i=1}^{n}{\bar{\mathcal{V}}_{i}}\left\\{K_{h}(S_{i}-s)Y_{i}\widehat{\omega}^{*}_{ai}-(\widehat{\omega}^{*}_{ai}-1)\widehat{\psi}^{*}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{X}_{i})\widehat{\psi}^{*}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{X}_{i})\right\\},$
$\displaystyle\widehat{f}^{*}_{\scriptscriptstyle a,\sf DR}(s)=$
$\displaystyle
n^{-1}\sum_{i=1}^{n}{\bar{\mathcal{V}}_{i}}\left\\{K_{h}(S_{i}-s)\widehat{\omega}^{*}_{ai}-(\widehat{\omega}^{*}_{ai}-1){\hat{\psi}_{a,f}^{*}(\mathbf{X}_{i})}\right\\},$

where

$\widehat{\psi}^{*}_{a,\scriptscriptstyle\sf
f}(s,\mathbf{x})=\frac{\sum_{i=1}^{n}\bar{\mathcal{V}}_{i}K_{\zeta}(\mathbf{X}_{i}^{\scriptscriptstyle\sf
T}\widehat{\boldsymbol{\gamma}}^{*}_{a}-\mathbf{x}^{\scriptscriptstyle\sf
T}\widehat{\boldsymbol{\gamma}}_{a}^{*})K_{h}(S_{i}-s)}{\sum_{i=1}^{n}\bar{\mathcal{V}}_{i}K_{\zeta}(\widehat{\boldsymbol{\gamma}}^{*}_{a}\mathbf{X}_{i}-\mathbf{x}^{\scriptscriptstyle\sf
T}\widehat{\boldsymbol{\gamma}}_{a}^{*})},\quad\widehat{\psi}^{*}_{a,\scriptscriptstyle\sf
m}(s,\mathbf{x})=M\\{\widehat{\boldsymbol{\beta}}_{a}^{*}(s)^{\scriptscriptstyle\sf
T}\vec{\mathbf{x}}\\}$

$\widehat{\boldsymbol{\gamma}}_{a}^{*}=\mbox{argmax}_{\boldsymbol{\gamma}\in\Omega}\\{\sum_{i\neq
j,A_{i}=A_{j}=a}\bar{\mathcal{V}}_{i}\bar{\mathcal{V}}_{j}I(\mathbf{X}_{i}^{\scriptscriptstyle\sf
T}\boldsymbol{\gamma}>\mathbf{X}_{j}^{\scriptscriptstyle\sf
T}\boldsymbol{\gamma})I(S_{i}>S_{j})\\}$ and
$\widehat{\boldsymbol{\beta}}_{a}^{*}(s)$ is the solution to

$\widehat{\mathbf{U}}_{a}^{*}(\boldsymbol{\beta};s)\equiv
n^{-1}\sum_{i=1}^{n}\bar{\mathcal{V}}_{i}I(A_{i}=a)K_{h}(S_{i}-s)\vec{\mathbf{X}}_{i}\left\\{Y_{i}-M(\boldsymbol{\beta}^{\scriptscriptstyle\sf
T}\vec{\mathbf{X}}_{i})\right\\}=\mathbf{0}.$

We construct the perturbed counterparts of $\widehat{g}_{\scriptscriptstyle\sf
DR}(s)$, $\hat{\Delta}_{\scriptscriptstyle\sf DR}$,
$\hat{\Delta}_{\widehat{g},\scriptscriptstyle\sf DR}$, and
$\widehat{\mbox{PTE}}_{\widehat{g},\scriptscriptstyle\sf DR}$ respectively as:

$\widehat{g}_{\scriptscriptstyle\sf
DR}^{*}(s)=\widehat{m}_{\scriptscriptstyle\sf
DR}^{*}(s)+\widehat{\lambda}_{\scriptscriptstyle\sf
DR}^{*}\widehat{\mathcal{P}}_{\scriptscriptstyle 0,\sf
DR}^{*}(s),\quad\widehat{\Delta}_{\scriptscriptstyle\sf
DR}^{*}=\widehat{\mu}_{\scriptscriptstyle 1,\sf
DR}^{*}-\widehat{\mu}_{\scriptscriptstyle 0,\sf
DR}^{*},\quad\widehat{\Delta}_{\widehat{g}^{*}_{\scriptscriptstyle\sf
DR},\scriptscriptstyle\sf
DR}^{*}=\widehat{\mu}_{1,\widehat{g}_{\scriptscriptstyle\sf
DR}^{*},\scriptscriptstyle\sf
DR}^{*}-\widehat{\mu}_{0,\widehat{g}_{\scriptscriptstyle\sf
DR}^{*},\scriptscriptstyle\sf DR}^{*},$

and $\widehat{\mbox{PTE}}_{\widehat{g}^{*}_{\scriptscriptstyle\sf
DR},\scriptscriptstyle\sf
DR}^{*}=\widehat{\Delta}_{\widehat{g}^{*}_{\scriptscriptstyle\sf
DR},\scriptscriptstyle\sf DR}^{*}/\widehat{\Delta}_{\scriptscriptstyle\sf
DR}^{*}$, where $\widehat{m}_{\scriptscriptstyle\sf
DR}^{*}(s)=\sum_{a=0}^{1}\widehat{m}_{\scriptscriptstyle a,\sf
DR}^{*}(s)\widehat{\mathcal{P}}^{*}_{\scriptscriptstyle a,\sf DR}(s)$,

$\widehat{\lambda}_{\scriptscriptstyle\sf
DR}^{*}=\frac{\int\left\\{\widehat{m}_{\scriptscriptstyle 0,\sf
DR}^{*}(s)-\widehat{m}^{*}_{\scriptscriptstyle 1,\sf
DR}(s)\right\\}\widehat{\mathcal{P}}_{\scriptscriptstyle 1,\sf
DR}^{*}(s)\widehat{f}^{*}_{\scriptscriptstyle 0,\sf
DR}(s)ds}{\int\widehat{\mathcal{P}}^{*}_{\scriptscriptstyle 0,\sf
DR}(s)\widehat{f}^{*}_{\scriptscriptstyle 0,\sf
DR}(s)ds},\quad\widehat{\mathcal{P}}^{*}_{\scriptscriptstyle a,\sf
DR}(s)=\frac{\widehat{f}^{*}_{\scriptscriptstyle a,\sf
DR}(s)}{\widehat{f}^{*}_{\scriptscriptstyle 0,\sf
DR}(s)+\widehat{f}^{*}_{\scriptscriptstyle 1,\sf DR}(s)},$
$\widehat{\mu}_{a,\widehat{g}_{\scriptscriptstyle\sf
DR}^{*},\scriptscriptstyle\sf
DR}^{*}=n^{-1}\sum_{i=1}^{n}\bar{\mathcal{V}}_{i}\left\\{\widehat{g}^{*}_{\scriptscriptstyle\sf
DR}(S_{i})\widehat{\omega}_{ai}^{*}-(\widehat{\omega}_{ai}^{*}-1)\hat{\zeta}_{a,\widehat{g}^{*}_{\scriptscriptstyle\sf
DR}}^{*}(\mathbf{X}_{i})\right\\},$ $\widehat{\mu}_{\scriptscriptstyle a,\sf
DR}^{*}=n^{-1}\sum_{i=1}^{n}\bar{\mathcal{V}}_{i}\left\\{Y_{i}\widehat{\omega}_{ai}^{*}-(\widehat{\omega}_{ai}^{*}-1)\hat{\zeta}_{a}^{*}(\mathbf{X}_{i})\right\\},$

$\hat{\zeta}_{a,g}^{*}(\mathbf{x})=\int
g(s)\widehat{\psi}_{a,\scriptscriptstyle\sf f}^{*}(s,\mathbf{x})ds$, and
$\hat{\zeta}_{a}^{*}(\mathbf{x})=\int\widehat{\psi}^{*}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{x})\widehat{\psi}^{*}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{x})ds$.

Operationally, we generate a large number, say $B=500$, realizations for
$\mathbf{V}$ and then obtain $B$ realizations of the perturbed statistics of
interest. Standard error estimates and confidence intervals can then be
constructed based on empirical variances of these realizations.

##  5 Simulation Studies

We have conducted simulation studies to evaluate the finite sample performance
of our proposed estimators compared to several existing methods. Namely, we
considered the naive PTE estimator of Freedman et al. (1992), denoted
$\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,naive}$, which does not take
into account baseline covariates $\mathbf{X}$ in the models $Y\mid S,A$ and
$Y\mid A$; a modified version that incorporates $\mathbf{X}$ into both models
$Y\mid S,A,\mathbf{X}$ and $Y\mid A,\mathbf{X}$, denoted
$\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,\mathbf{X}}$; a modified
version that incorporates the propensity score into both models, denoted
$\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,IPW}$; the PTE estimator given
in Parast et al. (2016), denoted $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf
P}$; and the PTE estimator of Wang et al. (2020), denoted
$\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf W}$. We let $n=400$ and $1000$
and choose $K(\cdot)$ as a Gaussian kernel. To obtain the bandwidth $h$ that
satisfies the undersmoothing assumption, we set $h=h_{opt}n^{-c_{0}}$,
$c_{0}=0.11$ where $h_{opt}=1.06n^{-1/5}$ as in Scott (2015). We compute the
true population parameters via Monte Carlo, under the counterfactual models
used to generate the data, with $N=100,000$ averaged over $100$ replications.
All results are summarized based on $500$ simulated datasets for each
configuration, and $B=500$ resampling replications were used for variance and
interval estimation based on the empirical variances.

We consider two general settings with a moderately strong surrogate in the
first setting and a weak surrogate in the second setting. Specifically, in the
first setting, we generate a 3-dimensional baseline covariate vector
$\mathbf{X}_{i}=(X_{i1},X_{i2},X_{i3})^{\scriptscriptstyle\sf T}$ as
$X_{i1}\sim\mathcal{N}(0,0.04)$, $X_{i2}\sim\mbox{Gamma}(2,2)$ and
$X_{i3}\sim\mbox{Uniform}(-1,1)$, and

$\displaystyle S_{i}^{(0)}$
$\displaystyle=\boldsymbol{\gamma}_{0[1]}^{\scriptscriptstyle\sf
T}\vec{\mathbf{X}}_{i}+{\epsilon_{i}},\quad
S_{i}^{(1)}=\boldsymbol{\gamma}_{1[1]}^{\scriptscriptstyle\sf
T}\vec{\mathbf{X}}_{i}+{\epsilon_{i}},$ $\displaystyle Y_{i}^{(0)}$
$\displaystyle=0.5S_{i}^{(0)}+\boldsymbol{\beta}_{0[1]}^{\scriptscriptstyle\sf
T}\vec{\mathbf{X}}_{i}+X_{1i}X_{2i}+X_{2i}X_{3i}+{e_{i}},$ $\displaystyle
Y_{i}^{(1)}$
$\displaystyle=0.3S_{i}^{(1)}+\boldsymbol{\beta}_{1[1]}^{\scriptscriptstyle\sf
T}\vec{\mathbf{X}}_{i}+X_{1i}X_{2i}+X_{2i}X_{3i}+{e_{i}},$

where $\epsilon_{i}\sim\mathcal{N}(0,1)$, $e_{i}\sim\mathcal{N}(0,0.04)$,
$\boldsymbol{\gamma}_{0[1]}=(0,0.5,1,-0.5)^{\scriptscriptstyle\sf T}$,
$\boldsymbol{\gamma}_{1[1]}=(0,1,0.5,2)^{\scriptscriptstyle\sf T}$,
$\boldsymbol{\beta}_{0[1]}=(0,0.2,-0.3,-0.5)^{\scriptscriptstyle\sf T}$, and
$\boldsymbol{\beta}_{1[1]}=(0,1,-0.5,0.2)^{\scriptscriptstyle\sf T}$. We
generate $A_{i}\mid\mathbf{X}_{i}$ from the propensity score model

$\displaystyle
P(A_{i}=1\mid\mathbf{X}_{i})=\mbox{expit}\\{-0.8X_{i1}+0.7X_{i2}-\log(X_{i3})+0.6X_{i1}X_{i3}\\}.$
(6)

Under this setting, $\Delta_{g_{\scriptscriptstyle\sf opt}}=0.29$ and
$\Delta=0.54$ so that the true potential outcomes PTE is $0.537$, i.e., $S$ is
a moderately strong surrogate for $Y$. We consider scenarios in which we
correctly specify both the PS and OR models, misspecify the PS model by
omitting the interaction term $X_{2}X_{3}$, misspecify the OR model by
omitting the variable $X_{2}$ and all interaction terms including $X_{2}$, and
misspecify both models.

In the second setting, we consider a relatively weak surrogate generate data
such that the effect of $S$ on $Y$ is non-linear. We generate baseline
covariates $\mathbf{X}_{i}=(X_{i1},X_{i2},X_{i3})^{\scriptscriptstyle\sf T}$
from $X_{i1}\sim\mathcal{N}(0,1)$, $X_{i2}\sim\mbox{Gamma}(2,2)$, and
$X_{i3}\sim\mbox{Uniform}(0,5)$. Given $\mathbf{X}_{i}$, we generate $S$ and
$Y$ from

$\displaystyle S_{i}^{(0)}=\boldsymbol{\gamma}_{0[2]}^{\scriptscriptstyle\sf
T}\vec{\mathbf{X}}_{i}+{\epsilon},$ $\displaystyle\quad
S_{i}^{(1)}=\boldsymbol{\gamma}_{1[2]}^{\scriptscriptstyle\sf
T}\vec{\mathbf{X}}_{i}+{\epsilon_{i}},$ $\displaystyle
Y_{i}^{(0)}=100+\boldsymbol{\beta}_{0[2]}(S_{i}^{(0)})^{\scriptscriptstyle\sf
T}\mathbf{X}_{i}+e_{i},$ $\displaystyle\quad
Y_{i}^{(1)}=50+\boldsymbol{\beta}_{1[2]}(S_{i}^{(1)})^{\scriptscriptstyle\sf
T}\mathbf{X}_{i}+e_{i},$

where $\epsilon_{i}\sim\mathcal{N}(0,4)$ and $e_{i}\sim\mathcal{N}(0,1)$ and
we let $\boldsymbol{\gamma}_{0[2]}=(100,1,5,0)^{\scriptscriptstyle\sf T}$,
$\boldsymbol{\gamma}_{1[2]}=(100,2,4,0)^{\scriptscriptstyle\sf T}$,
$\boldsymbol{\beta}_{0[2]}(s)=(s,-2\log(s),25)^{\scriptscriptstyle\sf T}$, and
$\boldsymbol{\beta}_{1[2]}(s)=(s,-3\log(s),-14)^{\scriptscriptstyle\sf T}$. We
generate $A_{i}\mid\mathbf{X}_{i}$ from model (6) as in setting 1. Under this
data generating mechanism, $\Delta_{g_{\scriptscriptstyle\sf opt}}=5.7$ and
$\Delta=26.7$, resulting in $\mbox{PTE}_{g_{\scriptscriptstyle\sf
opt}}=0.214$. We consider scenarios in which we correctly specify both the PS
and OR models, misspecify the PS model by omitting the $\log(X_{3})$ term,
misspecify the OR model by omitting $X_{2}$, and misspecify both models.

We first summarize results for setting 1. In Figure 1, we plot the empirical
biases, the empirical standard error (ESE) compared to the average of the
estimated standard error (ASE), and empirical coverage probabilities of the
95% pointwise confidence intervals (CIs) for $g_{opt}(\cdot)$ based on the DR
estimator $\widehat{g}_{\scriptscriptstyle\sf DR}(\cdot)$ estimated with
sample size $n=1000$ when (A) both the PS model and OR model are correctly
specified, (B) the PS model is misspecified but the OR model is correctly
specified, (C) the OR model is misspecified but the PS model is correctly
specified, and (D) both models are misspecified. When at least one of the two
models is correctly specified, the point estimates for $g_{opt}(\cdot)$
present negligible bias, the ASEs are close to the ESEs, probabilities of the
$95\%$ CIs are close to their nominal level. When both models are
misspecified, bias is observed in the tails, the ASE somewhat underestimates
the ESE, and the coverage probabilities of the $95\%$ confidence intervals are
somewhat below the nominal level. Results for $n=400$ bear similar patterns
and are hence omitted for brevity.

Figure 1: Empirical bias, empirical standard error (ESE) versus average of the
estimated standard error (ASE), and coverage probabilities of the 95%
confidence intervals for $\widehat{g}(s)$ when $n=1000$ and (Row 1) both
models are correctly specified, (Row 2) PS model is misspecified, (Row 3) OR
model is misspecified, (Row 4) both models are misspecified.

In Table 1, we summarize results for PTE estimation obtained via the proposed
method and other existing methods. When at least one of the PS and OR models
are correctly specified, the DR estimator displays negligible bias and nominal
coverage. The IPW estimator $\widehat{\mbox{PTE}}_{\widehat{g}}$ has
substantial bias when the PS model is misspecified and the DR estimator also
presents bias when both models are incorrect, as expected. In addition, the
IPW estimator is less efficient compared to the DR estimator when the PS model
is correctly specified. The estimate from
$\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf W}$, which assumes that treatment
is randomly assigned, shows considerable bias and below nominal coverage. It
is difficult to make a direct comparison of our estimators with other
literature estimators, as they estimate slightly different $PTE$ quantities,
but it is useful to note that $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf
F,\mathbf{X}}$ increases slightly from
$\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,naive}$ toward the true $PTE$
when baseline covariates are included correctly in the OR models $Y\mid
A,S,\mathbf{X}$ and $Y\mid A,\mathbf{X}$. However, when the OR models are
misspecified, $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,\mathbf{X}}$
shows substantial bias.

Table 1: Average of the estimated PTE, Bias for PTE estimators whose target parameter is $\mbox{PTE}_{g_{\scriptscriptstyle\sf opt}}=0.539$, Empirical Standard Error (ESE), Average of the Estimated Standard Errors (ASE), Root Mean Squared Error (RMSE), and Empirical Coverage Probabilities of the $95\%$ CIs of Estimators under Different Model Scenarios for setting 1. Size | Estimator | Scenario | Est PTE | Bias | ESE | ASE | RMSE | Coverage
---|---|---|---|---|---|---|---|---
n = 400 | $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,naive}$ |  No $\mathbf{X}$ | 0.380 | - | 0.114 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,\mathbf{X}}$ | OR Correct | 0.436 | - | 0.072 | - | - | -
|  | OR Misspecified | 0.611 | - | 0.122 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,IPW}$ | PS Correct | 0.460 | - | 0.203 | - | - | -
|  | PS Misspecified | 0.434 | - | 0.201 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf P}$ |  No $\mathbf{X}$ | 0.316 | - | 0.101 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf W}$ |  No $\mathbf{X}$ | 0.435 | -0.105 | 0.110 | 0.110 | 0.191 | 0.729
| $\widehat{\mbox{PTE}}_{\widehat{g}}$ | PS Correct | 0.545 | 0.006 | 0.087 | 0.088 | 0.088 | 0.930
|  | PS Misspecified | 0.480 | -0.059 | 0.107 | 0.109 | 0.122 | 0.916
| $\widehat{\mbox{PTE}}_{\widehat{g},\scriptscriptstyle\sf DR}$ | Both Correct | 0.542 | 0.003 | 0.079 | 0.079 | 0.080 | 0.940
|  | PS Misspecified | 0.534 | -0.005 | 0.074 | 0.079 | 0.075 | 0.954
|  | OR Misspecified | 0.532 | -0.007 | 0.084 | 0.082 | 0.085 | 0.940
|  | Both Misspecified | 0.448 | -0.091 | 0.119 | 0.115 | 0.149 | 0.779
n = 1000 | $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,naive}$ |  No $\mathbf{X}$ | 0.388 | - | 0.064 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,\mathbf{X}}$ | OR Correct | 0.445 | - | 0.040 | - | - | -
|  | OR Misspecified | 0.601 | - | 0.065 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,IPW}$ | PS Correct | 0.384 | - | 0.164 | - | - | -
|  | PS Misspecified | 0.372 | - | 0.159 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf P}$ |  No $\mathbf{X}$ | 0.335 | - | 0.061 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf W}$ |  No $\mathbf{X}$ | 0.437 | -0.103 | 0.181 | 0.174 | 0.240 | 0.803
| $\widehat{\mbox{PTE}}_{\widehat{g}}$ | PS Correct | 0.534 | -0.006 | 0.052 | 0.053 | 0.051 | 0.950
|  | PS Misspecified | 0.453 | -0.087 | 0.069 | 0.072 | 0.113 | 0.768
| $\widehat{\mbox{PTE}}_{\widehat{g},\scriptscriptstyle\sf DR}$ | Both Correct | 0.533 | -0.006 | 0.051 | 0.050 | 0.052 | 0.944
|  | PS Misspecified | 0.536 | -0.003 | 0.050 | 0.050 | 0.050 | 0.956
|  | OR Misspecified | 0.540 | 0.001 | 0.050 | 0.052 | 0.051 | 0.948
|  | Both Misspecified | 0.432 | -0.107 | 0.074 | 0.072 | 0.114 | 0.635

Table 2 shows that under setting 2, $\widehat{\mbox{PTE}}_{\widehat{g}}$ is
consistent when the PS model is correctly specified and
$\widehat{\mbox{PTE}}_{\widehat{g},\scriptscriptstyle\sf DR}$ is consistent
when either the PS model or OR model is correctly specified. However,
$\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,naive}$,
$\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,\mathbf{X}}$,
$\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,IPW}$, and
$\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf P}$ all estimate the true PTE as
being close to $0$. This is in part due to the nonmonotone relationship
between $Y$ and $S$ and the fact that these estimators use $S$ directly rather
than $g(S)$ in estimating the treatment effect.

Table 2: Estimated PTE, Bias for PTE estimators whose target parameter is $\mbox{PTE}_{g_{\scriptscriptstyle\sf opt}}=0.214$, Empirical Standard Error (ESE), Average of the Estimated Standard Errors (ASE), Root Mean Squared Error (RMSE),and and Empirical Coverage Probabilities of the $95\%$ CIs of Estimators under Different Model Scenarios for setting 2. Size | Estimator | Scenario | Est PTE | Bias | ESE | ASE | RMSE | Coverage
---|---|---|---|---|---|---|---|---
n = 400 | $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,naive}$ |  No $\mathbf{X}$ | 0.048 | - | 0.031 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,\mathbf{X}}$ | OR Correct | 0.063 | - | 0.051 | - | - | -
|  | OR Misspecified | 0.062 | - | 0.052 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,IPW}$ | PS Correct | 0.142 | - | 0.188 | - | - | -
|  | PS Misspecified | 0.155 | - | 0.220 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf P}$ |  No $\mathbf{X}$ | 0.035 | - | 0.044 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf W}$ |  No $\mathbf{X}$ | 0.199 | -0.015 | 0.042 | 0.041 | 0.042 | 0.980
| $\widehat{\mbox{PTE}}_{\widehat{g}}$ | PS Correct | 0.220 | 0.006 | 0.050 | 0.048 | 0.050 | 0.936
|  | PS Misspecified | 0.237 | 0.024 | 0.053 | 0.052 | 0.056 | 0.894
| $\widehat{\mbox{PTE}}_{\widehat{g},\scriptscriptstyle\sf DR}$ | Both Correct | 0.216 | 0.002 | 0.048 | 0.049 | 0.049 | 0.946
|  | PS Misspecified | 0.214 | 0.000 | 0.047 | 0.046 | 0.048 | 0.952
|  | OR Misspecified | 0.219 | 0.005 | 0.048 | 0.046 | 0.048 | 0.940
|  | Both Misspecified | 0.197 | -0.017 | 0.058 | 0.060 | 0.061 | 0.845
n = 1000 | $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,naive}$ |  No $\mathbf{X}$ | 0.004 | - | 0.005 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,\mathbf{X}}$ | OR Correct | 0.029 | - | 0.009 | - | - | -
|  | OR Misspecified | 0.028 | - | 0.009 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf F,IPW}$ | PS Correct | 0.153 | - | 0.182 | - | - | -
|  | PS Misspecified | 0.141 | - | 0.157 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf P}$ |  No $\mathbf{X}$ | -0.016 | - | 0.104 | - | - | -
| $\widehat{\mbox{PTE}}_{\scriptscriptstyle\sf W}$ |  No $\mathbf{X}$ | 0.202 | -0.011 | 0.032 | 0.030 | 0.031 | 0.970
| $\widehat{\mbox{PTE}}_{\widehat{g}}$ | PS Correct | 0.210 | -0.004 | 0.030 | 0.030 | 0.031 | 0.948
|  | PS Misspecified | 0.232 | 0.020 | 0.030 | 0.031 | 0.031 | 0.914
| $\widehat{\mbox{PTE}}_{\widehat{g},\scriptscriptstyle\sf DR}$ | Both Correct | 0.218 | 0.005 | 0.028 | 0.029 | 0.029 | 0.942
|  | PS Misspecified | 0.217 | 0.004 | 0.029 | 0.029 | 0.030 | 0.946
|  | OR Misspecified | 0.220 | 0.007 | 0.029 | 0.028 | 0.028 | 0.940
|  | Both Misspecified | 0.189 | -0.024 | 0.051 | 0.045 | 0.057 | 0.853

##  6 Data Application

Published randomized trials have shown that the partial Mayo score may be a
good surrogate for the full Mayo score in assessing biologic therapies for
ulcerative colitis (UC) (Lewis et al.,2008; Colombel et al., 2011;
Ananthakrishnan et al., 2016). The partial Mayo score is an inexpensive, non-
invasive composite score that can be measured early and ranges from $0$ to
$9$. It is based on a patient’s self-assessed stool frequency (0-3), rectal
bleeding (0-3), and a physician’s global assessment (0-3). The full Mayo score
ranges from $0$ to $12$ and consists of the partial Mayo score, in addition to
an invasive endoscopoy score evaluating mucosal appearance (0-3), and is
typically collected later in the trial. When conventional treatments fail,
biologic therapies such as infliximab, adalimumab, or golimumab may be used
(Rubin et al., 2019). While these medications have been shown to be effective
in placebo-controlled trials for rheumatoid arthritis (Taylor et al., 2017),
there has been a lack of trials comparing agents directly in UC. One of the
first such trials, a phase 3b trial of vedolizumab vs. adalimumab for patients
with moderate-to-severe UC, showed that vedolizumab was superior to adalimumab
with respect to achievement of clinical remission and endoscopic improvement,
but not corticosteroid-free clinical remission (Sands et al., 2019). The
researchers were unable to postulate an explanation for the inconsistency of
the results between the primary and secondary remission outcomes, and
concluded that this question required further investigation.

To illustrate the utility of our proposed methods, we apply our procedure to
examine the surrogacy of the partial Mayo score at week 6 on the primary
outcome of the full Mayo score at week 54 among patients with moderate-to-
severe UC. We examine an application of real-world interest in comparing head-
to-head trials of two biologic therapies for patients with UC (Ungaro &
Colombel, 2017). Treatment randomization is broken by combining data from two
separate trials on patients with active UC, one comparing infliximab against a
placebo (NCT00036439) and another comparing golimumab against a placebo
(NCT00488631). To adjust for confounding bias, we consider baseline covariates
$\mathbf{X}$ including patient age, sex, race, and a health status score
ranging from 0 to 100. Data was obtained from the Yale University Open Data
Access (YODA) database (Ross et al., 2018). The ranges of $S$ in the two
treatment groups are $\\{1,2,...,9\\}$, although the distributions are
somewhat different, as evidenced in Figure 3 (see Supplementary Materials).
The distributions of the primary outcome in the two treatment groups is
provided in Figure 4 (Supplementary Materials). The analysis focused on the
$381$ patients who had complete information on the partial Mayo score at week
6, the full Mayo score at week 54, and baseline covariates, with $216$
patients in the golimumab group and $165$ in the infliximab group. We applied
the proposed methods to examine $g_{opt}(\cdot)$ of the surrogate for
predicting the treatment response as quantified by the full Mayo score. The
estimated $g_{opt}(\cdot)$ along with point-wise CIs based on the IPW (red)
and DR (black) estimators are very similar. The estimated transformation
function appears to be slightly non-linear, although there is clearly a
positive trend between $s$ and $\widehat{g}_{opt}(s)$, as shown in Figure 2.

Figure 2: Estimated $g(s)$ based on IPW (red) and DR (black) estimators and
pointwise 95% confidence intervals for the partial Mayo score at week 6
(surrogate) in a cross-trial comparison of infliximab and golimumab for $361$
UC patients

The DR estimator for the treatment effect is estimated as
$\widehat{\Delta}=2.326$ in favor of golimumab and the corresponding treatment
effect on the predicted outcome $\Delta_{g_{\scriptscriptstyle\sf opt}}$ is
estimated to be $\widehat{\Delta}_{\widehat{g},\scriptscriptstyle\sf
DR}=2.031$. This results in a DR PTE estimate of $0.872$ with a $95\%$ CI of
$(0.741,1.003)$, suggesting that the partial Mayo score at week 6 is a strong
surrogate for the full Mayo score at week 54. The results for the IPW PTE
estimate is very similar at $0.878$ with a $95\%$ CI of $(0.728,1.029)$.

##  7 Discussion

There is great interest in leveraging RWD, including EHRs, registry data, and
cross-trial data, to inform the design of shorter and cheaper clinical trials
through surrogate marker validation. We propose the first IPW and DR
estimators for the PTE explained by a surrogate marker when treatment is not
randomly assigned. We generalize the approach detailed in Wang et al. (2020)
for RCT data to RWD in the presence of treatment by indication bias. Our
proposed doubly robust estimator is efficient and consistent when at least one
of the PS and OR models is correctly specified. In the case of UC, we have
validated a partial Mayo score at week 6, which does not require an invasive
endoscopy procedure, as a strong surrogate for the full Mayo score at week 54
in a cross-trial study, supplementing evidence from previous placebo-
controlled trials (Lewis et al., 2008; Colombel et al., 2011; Ananthakrishnan
et al., 2016). This finding may be particularly useful in informing future
cross-trial designs for biologic therapies.

To provide flexibility in the estimation of $g_{opt}(\cdot)$ and $PTE$, we use
a varying-coefficient model to estimate the conditional mean of
$Y_{i}^{(a)}\mid S_{i}^{(a)},\mathbf{X}_{i}$. We are able to handle multiple
confounders by implementing a two-step estimator that first reduces
potentially high-dimensional $\mathbf{X}$ into
$\mathbf{X}^{\scriptscriptstyle\sf T}\boldsymbol{\hat{\gamma}_{a}}$ through
the generalized regression model and then estimates the conditional density of
$S\mid\mathbf{X}^{\scriptscriptstyle\sf T}\boldsymbol{\hat{\gamma}_{a}}$ using
the method of Hall et al. (2004). This procedure may be computationally
intensive, and in the case when the researcher is confident in the
specification of the PS model, it may be advisable to consider the IPW
estimation procedure.

Our approach has some limitations. First, our proposed plug-in estimators for
PTE use the same data to estimate both $g_{opt}$ and PTE given $g$, which may
result in overfitting bias. However, in simulation studies, the bias appears
small compared to the standard error, even with modest sample sizes. For small
sample sizes, cross-validation may be needed, in which separate data is used
to estimate $g_{opt}$ and PTE given $g$. Second, our approach relies on a few
assumptions, the strongest of which is the working independence assumption
$(Y^{(1)},S^{(1)})\perp(Y^{(0)},S^{(0)})$ needed for deriving the form of
$g_{\scriptscriptstyle\sf opt}$. This assumption has been discussed
extensively in Wang et al. (2020). Here, we reiterate that the assumption is
only a working assumption that allows for derivation of the specific form of
$g_{opt}(\cdot)$ and is not required for valid inference. When the working
independence assumption is severely violated, our proposed $g_{opt}$ can still
be considered an optimal transformation of the surrogate marker for the
difference in the primary outcome for two independent patients, one in the
treatment group and the other in the control group. Third, we fit our PS
models using logistic regression with specified basis functions, but
alternative approaches like gradient boosting, super learner, and other
machine learning classifiers may be considered in future research (Parast &
Griffin, 2017).

## Acknowledgements

This study, carried out under YODA Project $\\#$ 2019-4092, used data obtained
from the Yale University Open Data Access Project, which has an agreement with
JANSSEN RESEARCH $\&$ DEVELOPMENT, L.L.C.. The interpretation and reporting of
research using this data are solely the responsibility of the authors and does
not necessarily represent the official views of the Yale University Open Data
Access Project or JANSSEN RESEARCH $\&$ DEVELOPMENT, L.L.C.. Larry Han was
supported by the Harvard Big Data Training Grant (NIH-funded T32) and the
Clinical Orthopedic and Musculoskeletal Education and Training (COMET) Program
(NIH-funded T32) housed at Brigham and Women’s Hospital, Harvard Medical
School, and Harvard T.H. Chan School of Public Health.

## Supplementary Materials

Supplementary material available online includes supplementary figures,
additional simulation results, and R code.

Figure 3: Histogram of the full Mayo score at week 54 (primary outcome) in the
two treatment groups Figure 4: Histogram of the partial Mayo score at week 6
(surrogate) in the two treatment groups

## References

  * (1) Ananthakrishnan, A.N., Cagan, A., Cai, T., Gainer, V.S., Shaw, S.Y., Sanova, G., Churchill, S., Karlson, E.W. Kohane, I., Liao, K.P. et al. (2016). Comparative effectiveness of infliximab and adalimumab in crohn’s disease and ulcerative colitis. Inflammatory bowel diseases 22, 880-885.
  * (2) Bentley, C., Cressman, S., van der Hoek, K., Arts, K., Dancey, J., & Peacock, S. (2019). Conducting clinical trials—costs, impacts, and the value of clinical trials networks: a scoping review. Clinical Trials, 16(2), 183-193.
  * (3) Ciani, O., Buyse, M., Drummond, M., Rasi, G., Saad, E. D., and Taylor, R. S. (2017). Time to review the role of surrogate end points in health policy: state of the art and the way forward. Value in Health, 20(3), 487-495.
  * (4) Colombel, J.F., Rutgeerts, P., Reinisch, W., Esser, D., Wang, Y., Lang, Y., Marano, C.W., Strauss, R., Oddens, B.J., Feagan, B.G. and Hanauer, S.B. (2011). Early mucosal healing with infliximab is associated with improved long-term clinical outcomes in ulcerative colitis. Gastroenterology, 141(4), 1194-1201
  * (5) Concato, J., Shah, N., and Horwitz, R. I. (2000). Randomized, controlled trials, observational studies, and the hierarchy of research designs. New England Journal of Medicine, 342(25), 1887-1892.
  * (6) Corrigan-Curay, J., Sacks, L., and Woodcock, J. (2018). Real-world evidence and real-world data for evaluating drug safety and effectiveness. Jama, 320(9), 867-868.
  * (7) Freedman, L. S., Graubard, B. I., and Schatzkin, A. (1992). Statistical validation of intermediate endpoints for chronic diseases. Statistics in medicine, 11(2), 167-178.
  * (8) Freedman, L. S., and Schatzkin, A. (1992). Sample size for studying intermediate endpoints within intervention trials or observational studies. American Journal of Epidemiology, 136(9), 1148-1159.
  * (9) Gyawali, B., Hey, S. P., and Kesselheim, A. S. (2020). Evaluating the evidence behind the surrogate measures included in the FDA’s table of surrogate endpoints as supporting approval of cancer drugs. EClinicalMedicine, 21, 100332.
  * (10) Hall, P., Racine, J., and Li, Q. (2004). Cross-validation and the estimation of conditional probability densities. Journal of the American Statistical Association, 99(468), 1015-1026.
  * (11) Han, A. K. (1987). Non-parametric analysis of a generalized regression model: the maximum rank correlation estimator. Journal of Econometrics, 35(2-3), 303-316.
  * (12) Hastie, T., and Tibshirani, R. (1993). Varying‐coefficient models. Journal of the Royal Statistical Society: Series B (Methodological), 55(4), 757-779.
  * (13) Hernán, M. A., and Robins, J. M. (2016). Using big data to emulate a target trial when a randomized trial is not available. American journal of epidemiology, 183(8), 758-764.
  * (14) Hey, S. P., Kesselheim, A. S., Patel, P., Mehrotra, P., and Powers, J. H. (2020). US Food and Drug Administration recommendations on the use of surrogate measures as end points in new anti-infective drug approvals. JAMA internal medicine, 180(1), 131-138.
  * (15) Imbens, G. W., and Rubin, D. B. (2015). Causal inference in statistics, social, and biomedical sciences. Cambridge University Press.
  * (16) Jin, Z., Ying, Z., and Wei, L. J. (2001). A simple resampling method by perturbing the minimand. Biometrika, 88(2), 381-390.
  * (17) Lewis, J. D., Chuai, S., Nessel, L., Lichtenstein, G. R., Aberra, F. N., and Ellenberg, J. H. (2008). Use of the noninvasive components of the Mayo score to assess clinical response in ulcerative colitis. Inflammatory bowel diseases, 14(12), 1660-1666.
  * (18) Lin, D. Y., Fleming, T. R., and De Gruttola, V. (1997). Estimating the proportion of treatment effect explained by a surrogate marker. Statistics in medicine, 16(13), 1515-1527.
  * (19) Lurie, N., Saville, M., Hatchett, R., and Halton, J. (2020). Developing Covid-19 vaccines at pandemic speed. New England Journal of Medicine, 382(21), 1969-1973.
  * (20) Masry, E. (1996). Multivariate local polynomial regression for time series: uniform strong consistency and rates. Journal of Time Series Analysis, 17(6), 571-599.
  * (21) Pagan, A., and Ullah, A. (1999). Nonparametric econometrics. Cambridge University Press.
  * (22) Parast, L., and Griffin, B. A. (2017). Landmark estimation of survival and treatment effects in observational studies. Lifetime data analysis, 23(2), 161-182.
  * (23) Parast, L., McDermott, M. M., and Tian, L. (2016). Robust estimation of the proportion of treatment effect explained by surrogate marker information. Statistics in medicine, 35(10), 1637-1653.
  * (24) rentice, R. L. (1989). Surrogate endpoints in clinical trials: definition and operational criteria. Statistics in medicine, 8(4), 431-440.
  * (25) Price, B. L., Gilbert, P. B., and van der Laan, M. J. (2018). Estimation of the optimal surrogate based on a randomized trial. Biometrics, 74(4), 1271-1281.
  * (26) Robins, J. M., Rotnitzky, A., and Zhao, L. P. (1994). Estimation of regression coefficients when some regressors are not always observed. Journal of the American statistical Association, 89(427), 846-866.
  * (27) Ross, J.S., Waldstreicher, J., Bamford, S., Berlin, J.A., Childers, K., Desai, N.R., Gamble, G., Gross, C.P., Kuntz, R., Lehman, R. and Lins, P. (2018). Overview and experience of the YODA Project with clinical trial data sharing after 5 years. Scientific data, 5(1), 1-14.
  * (28) Rubin, D. B. (2005). Causal inference using potential outcomes: Design, modeling, decisions. Journal of the American Statistical Association, 100(469), 322-331.
  * (29) Rubin, D. T., Ananthakrishnan, A. N., Siegel, C. A., Sauer, B. G., and Long, M. D. (2019). ACG clinical guideline: ulcerative colitis in adults. American Journal of Gastroenterology, 114(3), 384-413.
  * (30) Sands, B.E., Peyrin-Biroulet, L., Loftus Jr, E.V., Danese, S., Colombel, J.F., Törüner, M., Jonaitis, L., Abhyankar, B., Chen, J., Rogers, R. and Lirio, R.A. (2019). Vedolizumab versus adalimumab for moderate-to-severe ulcerative colitis. New England Journal of Medicine, 381(13), 1215-1226.
  * (31) Scott, D. W. (2015). Multivariate density estimation: theory, practice, and visualization. John Wiley & Sons.
  * (32) Sherman, R. P. (1993). The limiting distribution of the maximum rank correlation estimator. Econometrica: Journal of the Econometric Society, 123-137.
  * (33) Taylor, P.C., Keystone, E.C., van der Heijde, D., Weinblatt, M.E., del Carmen Morales, L., Reyes Gonzaga, J., Yakushin, S., Ishii, T., Emoto, K., Beattie, S. and Arora, V. (2017). Baricitinib versus placebo or adalimumab in rheumatoid arthritis. New England Journal of Medicine, 376(7), 652-662.
  * (34) Tian, L., Zucker, D., and Wei, L. J. (2005). On the Cox model with time-varying regression coefficients. Journal of the American statistical Association, 100(469), 172-183.
  * (35) Ungaro, R. C., and Colombel, J. F. (2017). Biologics in inflammatory bowel disease—time for direct comparisons. Alimentary pharmacology & therapeutics, 46(1), 68.
  * (36) VanderWeele, T. J. (2013). Surrogate measures and consistent surrogates. Biometrics, 69(3), 561-565.
  * (37) Wang, X., Parast, L., Tian, L., and Cai, T. (2020). Model-free approach to quantifying the proportion of treatment effect explained by a surrogate marker. Biometrika, 107(1), 107-122.
  * (38) Wang, Y., and Taylor, J. M. (2002). A measure of the proportion of treatment effect explained by a surrogate marker. Biometrics, 58(4), 803-812.
  * (39) Wickström, K., and Moseley, J. (2017). Biomarkers and surrogate endpoints in drug development: a European regulatory view. Investigative ophthalmology & visual science, 58(6), BIO27-BIO33.

## Appendix 1

### Consistency and asymptotic normality of
$\widehat{\mbox{PTE}}_{\widehat{g}}$

Throughout, we assume that all components of $(Y,S,A,\mathbf{X})$ are sub-
gaussian, the true conditional mean function $\psi_{a,\scriptscriptstyle\sf
m}^{\scriptscriptstyle{\dagger}}(s;\mathbf{x})=E(Y^{(a)}\mid
S^{(a)}=s,\mathbf{X}=\mathbf{x})$ and the true conditional density of
$S^{(a)}\mid\mathbf{X}=\mathbf{x}$, $\psi_{a,\scriptscriptstyle\sf
f}^{\scriptscriptstyle{\dagger}}(s;\mathbf{x})$, are continuously
differentiable. We also assume that $S^{(a)}$ has a finite support and that
$h=O(n^{-\nu})$ with $\nu\in(1/4,1/2)$. In this section, we show that when the
propensity score model is correctly specified, the proposed IPW kernel
smoothed estimators $\widehat{g}(s)$ and $\widehat{\mbox{PTE}}_{\widehat{g}}$
are consistent for $g_{\scriptscriptstyle\sf opt}(s)$ and
$\mbox{PTE}_{g_{\scriptscriptstyle\sf opt}}$, respectively. We will also show
that
$\sqrt{n}(\widehat{\mbox{PTE}}_{\widehat{g}}-\mbox{PTE}_{g_{\scriptscriptstyle\sf
opt}})$ converges in distribution to a normal distribution with mean zero and
variance $\sigma^{2}$, which we will derive. To this end, we first show that
$\widehat{m}_{a}(s)$ and $\widehat{f}_{a}(s)$ are consistent for $m_{a}(s)$
and $f_{a}(s)$, respectively. Without loss of generality, we prove the
consistency of
$\widehat{m}_{a}(s)\equiv\widehat{m}_{a}(s;\widehat{\boldsymbol{\alpha}})$ for
$m_{a}(s)=E(Y^{(a)}\mid S^{(a)}=s)=E\\{\psi_{a,\scriptscriptstyle\sf
m}^{\scriptscriptstyle{\dagger}}(s;\mathbf{X})\\}$, where

$\widehat{m}_{a}(s;\boldsymbol{\alpha})=\frac{n^{-1}\sum_{i=1}^{n}K_{h}(S_{i}-s)Y_{i}I(A_{i}=a)/\pi_{a}(\mathbf{X}_{i};\boldsymbol{\alpha})}{n^{-1}\sum_{i=1}^{n}K_{h}(S_{i}-s)I(A_{i}=a)/\pi_{a}(\mathbf{X}_{i};\boldsymbol{\alpha})}.$

First, under the correct specification of the PS model,
$\widehat{\boldsymbol{\alpha}}\to\boldsymbol{\alpha}_{0}$ in probability,
where $\boldsymbol{\alpha}_{0}$ is the true parameter value. Hence
$\max_{i}|\widehat{\omega}_{ai}-\omega_{ai}|\leq\sup_{\mathbf{x}}|\pi_{a}(\mathbf{x};\widehat{\boldsymbol{\alpha}})^{-1}-\pi_{a}(\mathbf{x};\boldsymbol{\alpha}_{0})^{-1}|\to
0$ in probability, where
$\omega_{ai}=I(A_{i}=a)/\pi_{a}(\mathbf{X}_{i};\boldsymbol{\alpha}_{0})$. It
then follows from standard theory for non-parametric kernel estimators (Masry,
1996; Pagan & Ullah, 1999) and Taylor series expansions that

$\displaystyle\sup_{s}|\widehat{m}_{a}(s;\widehat{\boldsymbol{\alpha}})-m_{a}(s)|\leq$
$\displaystyle\sup_{s}|\widehat{m}_{a}(s;\boldsymbol{\alpha}_{0})-m_{a}(s)|+\sup_{s,\boldsymbol{\alpha}:\|\boldsymbol{\alpha}-\boldsymbol{\alpha}_{0}\|\leq
c}\|\widehat{\mathbf{m}}^{\prime}_{a}(s;\boldsymbol{\alpha})\|_{2}\|\widehat{\boldsymbol{\alpha}}-\boldsymbol{\alpha}_{0}\|_{2}$
$\displaystyle=$ $\displaystyle O_{p}\\{(nh)^{-\frac{1}{2}}\sqrt{\log
n}+h^{2}+n^{-\frac{1}{2}}\\}=o_{p}(1)$

where
$\widehat{\mathbf{m}}_{a}^{\prime}(s;\boldsymbol{\alpha})=\partial\widehat{m}_{a}(s;\boldsymbol{\alpha})/\partial\boldsymbol{\alpha}$
and $c$ is any small constant. Similarly, we have
$\sup_{s}|\widehat{f}_{a}(s)-f_{a}(s)|=O_{p}\\{(nh)^{-\frac{1}{2}}\sqrt{\log
n}+h^{2}+n^{-\frac{1}{2}}\\}=o_{p}(1)$. When $h=O(n^{-\nu})$ with
$\nu\in(1/4,1/2)$, it is not difficult to show that
$\widehat{\lambda}-\lambda=O_{p}(n^{-\frac{1}{2}}+h^{2})=O_{p}(n^{-\frac{1}{2}})$.
It follows that
$\sup_{s}|\widehat{g}(s)-g(s)|=O_{p}\\{(nh)^{-\frac{1}{2}}\sqrt{\log
n}+h^{2}+n^{-\frac{1}{2}}\\}=o_{p}(1)$. Similarly, we may show that

$\displaystyle|\widehat{\Delta}_{\widehat{g}}-\Delta_{g_{\scriptscriptstyle\sf
opt}}|=\widehat{\Delta}_{\widehat{g}}-\widehat{\Delta}_{g_{\scriptscriptstyle\sf
opt}}+\widehat{\Delta}_{g_{\scriptscriptstyle\sf
opt}}-\Delta_{g_{\scriptscriptstyle\sf
opt}}=\int\\{\widehat{g}(s)-g_{\scriptscriptstyle\sf
opt}(s)\\}d\widehat{D}(s)+O_{p}(n^{-\frac{1}{2}}),$

where
$\widehat{D}(s)=n^{-1}\sum_{i=1}^{n}(\widehat{\omega}_{1i}-\widehat{\omega}_{0i})I(S_{i}\leq
s)$. It follows from the uniform convergence of $\widehat{g}(s)\to
g_{\scriptscriptstyle\sf opt}(s)$ and $\widehat{D}(s)\to D(s)=P(S^{(1)}\leq
s)-P(S^{(0)}\leq s)$ that
$\widehat{\Delta}_{\widehat{g}}-\Delta_{g_{\scriptscriptstyle\sf opt}}\to 0$
in probability. This, together with the consistency of $\widehat{\Delta}$ for
$\Delta$, implies the consistency of $\widehat{\mbox{PTE}}_{\widehat{g}}$ for
$\mbox{PTE}_{g_{\scriptscriptstyle\sf opt}}$.

We next establish the asymptotic normality of
$\sqrt{n}(\widehat{\mbox{PTE}}_{\widehat{g}}-\mbox{PTE}_{g_{\scriptscriptstyle\sf
opt}})$. First, note that

$\displaystyle\widehat{f}_{a}(s)-f_{a}(s)$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\left[\omega_{ai}\left\\{K_{h}\left(S_{i}-s\right)-f_{a}(s)\right\\}+\mathbf{f}^{\prime}_{a}(s)\boldsymbol{\mathcal{U}}_{\alpha,i}\right]+o_{p}((nh)^{-1/2}),$
$\displaystyle\widehat{m}_{a}(s)-m_{a}(s)$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\left[\omega_{ai}K_{h}\left(S_{i}-s\right)\mathcal{U}_{m_{a},i}(s)+\mathbf{m}^{\prime}_{a}(s)\boldsymbol{\mathcal{U}}_{\alpha,i}\right]+o_{p}((nh)^{-1/2}),$

where $\mathbf{m}_{a}^{\prime}(s;\boldsymbol{\alpha})=\partial
m_{a}(s;\boldsymbol{\alpha})/\partial\boldsymbol{\alpha}$,
$\mathbf{f}^{\prime}_{a}(s;\boldsymbol{\alpha})=\partial
f_{a}(s;\boldsymbol{\alpha})/\partial\boldsymbol{\alpha}$,
$\mathcal{U}_{m_{a},i}(s)=f_{a}(s)^{-1}\\{Y_{i}^{(a)}-m_{a}(s)\\}$,

$\displaystyle m_{a}(s;\boldsymbol{\alpha})$
$\displaystyle=\frac{E\left\\{\psi_{a,\scriptscriptstyle\sf
m}^{\scriptscriptstyle{\dagger}}(s;\mathbf{X}_{i})\psi_{a,\scriptscriptstyle\sf
f}^{\scriptscriptstyle{\dagger}}(s;\mathbf{X}_{i})\frac{\pi(\mathbf{X}_{i};\boldsymbol{\alpha}_{0})}{\pi(\mathbf{X}_{i};\boldsymbol{\alpha})}\right\\}}{f_{a}(s;\boldsymbol{\alpha})},\quad
f_{a}(s;\boldsymbol{\alpha})=E\left\\{\psi_{a,\scriptscriptstyle\sf
f}^{\scriptscriptstyle{\dagger}}(s;\mathbf{X}_{i})\frac{\pi(\mathbf{X}_{i};\boldsymbol{\alpha}_{0})}{\pi(\mathbf{X}_{i};\boldsymbol{\alpha})}\right\\},$

and
$\widehat{\boldsymbol{\alpha}}-\boldsymbol{\alpha}_{0}=n^{-1}\sum_{i=1}^{n}\boldsymbol{\mathcal{U}}_{\alpha,i}+o_{p}(n^{-\frac{1}{2}})$
following standard likelihood theory. It follows that

$\displaystyle\widehat{\mathcal{P}}_{0}(s)-\mathcal{P}_{0}(s)=\frac{\hat{f}_{0}(s)f_{1}(s)-\hat{f}_{1}(s)f_{0}(s)}{\left\\{f_{1}(s)+f_{0}(s)\right\\}^{2}}+o_{p}((nh)^{-1/2})$
$\displaystyle=\frac{\mathcal{P}_{1}(s)\mathcal{P}_{0}(s)(\hat{f}_{0}(s)-f_{0}(s))}{f_{0}(s)}-\frac{\mathcal{P}_{1}(s)\mathcal{P}_{0}(s)(\hat{f}_{1}(s)-f_{1}(s))}{f_{1}(s)}+o_{p}((nh)^{-1/2})$
$\displaystyle=\mathcal{P}_{1}(s)\mathcal{P}_{0}(s)n^{-1}\sum_{i=1}^{n}\left[\omega_{0i}\bigg{\\{}K_{h}(S_{i}-s)f_{0}(s)^{-1}-1\right\\}-\omega_{1i}\left\\{K_{h}(S_{i}-s)f_{1}(s)^{-1}-1\right\\}$
$\displaystyle+\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\left\\{\frac{\mathbf{f}^{\prime}_{0}(s)}{f_{0}(s)}-\frac{\mathbf{f}^{\prime}_{1}(s)}{f_{1}(s)}\right\\}\bigg{]}+o_{p}((nh)^{-1/2})$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\left[K_{h}(S_{i}-s)\mathcal{G}_{\mathcal{P}_{0},i}(s)+(\omega_{1i}-\omega_{0i})\mathcal{P}_{1}(s)\mathcal{P}_{0}(s)+\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\boldsymbol{\mathcal{B}}_{\mathcal{P}_{0}}(s)\right]+o_{p}((nh)^{-1/2})$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\mathcal{U}_{\mathcal{P}_{0},i}+o_{p}((nh)^{-1/2}),$

where

$\displaystyle\mathcal{U}_{\mathcal{P}_{0},i}(s)$
$\displaystyle=K_{h}(S_{i}-s)\mathcal{G}_{\mathcal{P}_{0},i}(s)+(\omega_{1i}-\omega_{0i})\mathcal{P}_{1}(s)\mathcal{P}_{0}(s)+\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\boldsymbol{\mathcal{B}}_{\mathcal{P}_{0}}(s),$
$\displaystyle\mathcal{G}_{\mathcal{P}_{0},i}(s)$
$\displaystyle=\\{\omega_{0i}/f_{0}(s)-\omega_{1i}/f_{1}(s)\\}\mathcal{P}_{1}(s)\mathcal{P}_{0}(s),$
$\displaystyle\boldsymbol{\mathcal{B}}_{\mathcal{P}_{0}}(s)$
$\displaystyle=\\{\mathbf{f}^{\prime}_{0}(s)/f_{0}(s)-\mathbf{f}^{\prime}_{1}(s)/f_{1}(s)\\}\mathcal{P}_{1}(s)\mathcal{P}_{0}(s)$

Similarly, we have
$\widehat{\mathcal{P}}_{1}(s)-\mathcal{P}_{1}(s)=-(\widehat{\mathcal{P}}_{0}(s)-\mathcal{P}_{0}(s))=-n^{-1}\sum_{i=1}^{n}\mathcal{U}_{\mathcal{P}_{0},i}(s)+o_{p}(n^{-\frac{1}{2}})$.

Now

$\displaystyle\widehat{m}(s)-m(s)$
$\displaystyle=\\{\widehat{m}_{1}(s)-m_{1}(s)\\}\mathcal{P}_{1}(s)+\widehat{m}_{1}(s)\\{\widehat{\mathcal{P}}_{1}(s)-\mathcal{P}_{1}(s)\\}$
$\displaystyle\quad+\\{\widehat{m}_{0}(s)-m_{0}(s)\\}\mathcal{P}_{0}(s)+\widehat{m}_{0}(s)\\{\widehat{\mathcal{P}}_{0}(s)-\mathcal{P}_{0}(s)\\}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\mathcal{U}_{m,i}(s)+o_{p}((nh)^{-1/2}),$

where

$\displaystyle\mathcal{U}_{m,i}(s)$
$\displaystyle=\sum_{a=0}^{1}\left\\{\omega_{ai}K_{h}(S_{i}-s)\mathcal{U}_{m_{a},i}(s)+\mathbf{m}_{a}^{\prime}(s)\boldsymbol{\mathcal{U}}_{\alpha,i}\right\\}\mathcal{P}_{a}(s)+\\{m_{0}(s)-m_{1}(s)\\}\mathcal{U}_{\mathcal{P}_{0},i}(s)$
$\displaystyle=K_{h}(S_{i}-s)\mathcal{G}_{m,i}(s)+\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\boldsymbol{\mathcal{B}}_{m}(s)+(\omega_{1i}-\omega_{0i})\mathcal{A}_{m}(s),$

where

$\displaystyle\mathcal{G}_{m,i}(s)$
$\displaystyle=\sum_{a=0}^{1}\omega_{ai}\mathcal{U}_{ma,i}(s)\mathcal{P}_{a}(s)+\mathcal{G}_{\mathcal{P}_{0},i}(s)\\{m_{0}(s)-m_{1}(s)\\},$
$\displaystyle\boldsymbol{\mathcal{B}}_{m}(s)$
$\displaystyle=\sum_{a=0}^{1}\mathbf{m}_{a}^{\prime}(s)\mathcal{P}_{a}(s)+\boldsymbol{\mathcal{B}}_{\mathcal{P}_{0}}(s)\\{m_{0}(s)-m_{1}(s)\\},$
$\displaystyle\mathcal{A}_{m}(s)$
$\displaystyle=\mathcal{P}_{1}(s)\mathcal{P}_{0}(s)\\{m_{0}(s)-m_{1}(s)\\}$

Together with arguments given in Appendix D of Wang et al. (2020), the fact
that $h=o_{p}(n^{-1/4})$, and a Taylor series expansion for approximating
$\int K_{h}(S_{i}-s)H(s)ds$ for any given smooth function $H$, where we denote
$\mu_{H_{0}}:=\int H(s)dF_{0}(s)$, we have the following expansion for
$\int\widehat{m}(s)d\widehat{F}_{0}(s)-\mu_{m0}$:

$\displaystyle\int\widehat{m}(s)\widehat{f}_{0}(s)ds-\mu_{m0}$
$\displaystyle=\int\\{\widehat{m}(s)-m(s)\\}f_{0}(s)ds$
$\displaystyle\quad+\int
m(s)\left\\{\widehat{f}_{0}(s)-f_{0}(s)\right\\}ds+o_{p}(n^{-\frac{1}{2}})$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\bigg{[}\int
K_{h}(S_{i}-s)\left\\{\mathcal{G}_{m,i}(s)f_{0}(s)+\omega_{0i}m(s)\right\\}ds$
$\displaystyle\quad+(\omega_{1i}-\omega_{0i})\int\mathcal{A}_{m}(s)f_{0}(s)ds-\omega_{0i}\int
m(s)f_{0}(s)ds$
$\displaystyle\quad+\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\int\left\\{m(s)\mathbf{f}^{\prime}_{0}(s)+\boldsymbol{\mathcal{B}}_{m}(s)f_{0}(s)\right\\}ds\bigg{]}+o_{p}(n^{-1/2})$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\bigg{[}\mathcal{G}_{m,i}(S_{i})f_{0}(S_{i})+\omega_{0i}m(S_{i})+(\omega_{1i}-\omega_{0i})\int\mathcal{A}_{m}(s)f_{0}(s)ds$
$\displaystyle\quad-\omega_{0i}\int
m(s)f_{0}(s)ds+\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\int\left\\{m(s)\mathbf{f}^{\prime}_{0}(s)+\boldsymbol{\mathcal{B}}_{m}(s)f_{0}(s)\right\\}ds\bigg{]}$
$\displaystyle\quad+o_{p}(n^{-1/2}).$

Similarly, we have

$\displaystyle\int\widehat{P}_{0}(s)$
$\displaystyle\widehat{f}_{0}(s)ds-\mu_{\mathcal{P}_{0}0}$
$\displaystyle=\int\left\\{\widehat{P}_{0}(s)-\mathcal{P}_{0}(s)\right\\}f_{0}(s)ds+\int\mathcal{P}_{0}(s)\left\\{\widehat{f}_{0}(s)-f_{0}(s)\right\\}ds+o_{p}(n^{-1/2})$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\bigg{[}\int
K_{h}(S_{i}-s)\\{\mathcal{G}_{\mathcal{P}_{0},i}(s)f_{0}(s)+\omega_{0i}\mathcal{P}_{0}(s)\\}ds$
$\displaystyle\quad+(\omega_{1i}-\omega_{0i})\int\mathcal{P}_{1}(s)\mathcal{P}_{0}(s)f_{0}(s)ds-\omega_{0i}\int\mathcal{P}_{0}(s)f_{0}(s)ds$
$\displaystyle\quad+\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\int\left\\{\mathcal{P}_{0}(s)\mathbf{f}^{\prime}_{0}(s)+\boldsymbol{\mathcal{B}}_{\mathcal{P}_{0}}(s)f_{0}(s)\right\\}ds\bigg{]}+o_{p}(n^{-1/2})$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\bigg{[}\mathcal{G}_{\mathcal{P}_{0},i}(S_{i})f_{0}(S_{i})+\omega_{0i}\mathcal{P}_{0}(S_{i})+(\omega_{1i}-\omega_{0i})\int\mathcal{P}_{1}(s)\mathcal{P}_{0}(s)f_{0}(s)ds$
$\displaystyle\quad-\omega_{0i}\int\mathcal{P}_{0}(s)f_{0}(s)ds+\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\int\left\\{\mathcal{P}_{0}(s)\mathbf{f}^{\prime}_{0}(s)+\boldsymbol{\mathcal{B}}_{\mathcal{P}_{0}}(s)f_{0}(s)\right\\}ds\bigg{]}$
$\displaystyle\quad+o_{p}(n^{-1/2}).$

Since
$\lambda=\frac{\displaystyle{\mu_{0}-\mu_{m0}}}{\displaystyle{\mu_{\mathcal{P}_{0}0}}}$
and
$\widehat{\lambda}=\frac{\displaystyle{\widehat{\mu}_{0}-\widehat{\mu}_{m0}}}{\displaystyle{\widehat{\mu}_{\mathcal{P}_{0}0}}}$,
it follows from above that

$\widehat{\lambda}-\lambda=n^{-1}\sum_{i=1}^{n}\mathcal{U}_{\lambda,i}+o_{p}(n^{-\frac{1}{2}}),$

where

$\displaystyle\widehat{\lambda}-\lambda$
$\displaystyle=\mu_{\mathcal{P}_{0}0}^{-1}\left(\widehat{\mu}_{0}-\mu_{0}\right)-\mu_{\mathcal{P}_{0}0}^{-1}\lambda\left\\{\int\widehat{\mathcal{P}}_{0}(s)d\widehat{F}_{0}(s)-\mu_{\mathcal{P}_{0}0}\right\\}$
$\displaystyle\quad-\mu_{\mathcal{P}_{0}0}^{-1}\left\\{\int\widehat{m}(s)d\widehat{F}_{0}(s)-\mu_{m0}\right\\}+o_{p}(n^{-1/2})$
$\displaystyle=\mu_{\mathcal{P}_{0}0}^{-1}n^{-1}\sum_{i=1}^{n}w_{0i}\left(Y_{i}-\mu_{0}\right)-\mu_{\mathcal{P}_{0}0}^{-1}\lambda
n^{-1}\sum_{i=1}^{n}\bigg{\\{}\mathcal{G}_{\mathcal{P}_{0},i}(S_{i})f_{0}(S_{i})+\omega_{0i}\mathcal{P}_{0}(S_{i})$
$\displaystyle\quad+(\omega_{1i}-\omega_{0i})\int\mathcal{P}_{1}(s)\mathcal{P}_{0}(s)f_{0}(s)ds-\omega_{0i}\int\mathcal{P}_{0}(s)f_{0}(s)ds\bigg{\\}}$
$\displaystyle\quad+\mu_{\mathcal{P}_{0}0}^{-1}\lambda
n^{-1}\sum_{i=1}^{n}\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\int\left\\{\mathcal{P}_{0}(s)\mathbf{f}^{\prime}_{0}(s)+\boldsymbol{\mathcal{B}}_{\mathcal{P}_{0}}(s)f_{0}(s)\right\\}ds$
$\displaystyle\quad-\mu_{\mathcal{P}_{0}0}^{-1}n^{-1}\sum_{i=1}^{n}\bigg{\\{}\mathcal{G}_{m,i}(S_{i})f_{0}(S_{i})+\omega_{0i}m(S_{i})$
$\displaystyle\quad+(\omega_{1i}-\omega_{0i})\int\mathcal{A}_{m}(s)f_{0}(s)ds-\omega_{0i}\int
m(s)f_{0}(s)ds\bigg{\\}}$
$\displaystyle\quad+\mu_{\mathcal{P}_{0}0}^{-1}n^{-1}\sum_{i=1}^{n}\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\int\left\\{m(s)\mathbf{f}^{\prime}_{0}(s)ds+\boldsymbol{\mathcal{B}}_{m}(s)f_{0}(s)\right\\}ds+o_{p}(n^{-1/2})$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\mathcal{U}_{\lambda}(\textbf{D}_{i})+o_{p}(n^{-1/2}).$

Gathering the above expansions, we may obtain the form of
$\widehat{g}(s)-g_{\scriptscriptstyle\sf opt}(s)$ as

$\displaystyle\widehat{g}(s)$ $\displaystyle-g_{\mathrm{opt}}(s)$
$\displaystyle=\widehat{m}(s)-m(s)+(\widehat{\lambda}-\lambda)\mathcal{P}_{0}(s)+\lambda\left\\{\widehat{\mathcal{P}}_{0}(s)-\mathcal{P}_{0}(s)\right\\}+o_{p}((nh)^{-1/2})$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\bigg{[}K_{h}(S_{i}-s)\mathcal{G}_{m,i}(s)+\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\boldsymbol{\mathcal{B}}_{m}(s)+(\omega_{1i}-\omega_{0i})\mathcal{A}_{m}(s)+\mathcal{P}_{0}(s)\mathcal{U}_{\lambda}(\textbf{D}_{i})$
$\displaystyle\quad+\lambda\\{K_{h}(S_{i}-s)\mathcal{G}_{\mathcal{P}_{0},i}(s)+(\omega_{1i}-\omega_{0i})\mathcal{P}_{1}(s)\mathcal{P}_{0}(s)+\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\boldsymbol{\mathcal{B}}_{\mathcal{P}_{0}}(s)\\}\bigg{]}$
$\displaystyle\quad+o_{p}((nh)^{-1/2})$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\mathcal{U}_{G}(s;\textbf{D}_{i})+o_{p}((nh)^{-1/2}),$

where

$\displaystyle\mathcal{U}_{G}(s;\textbf{D}_{i})$
$\displaystyle=K_{h}(S_{i}-s)\\{\mathcal{G}_{m,i}(s)+\lambda\mathcal{G}_{\mathcal{P}_{0},i}(s)\\}+\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\\{\boldsymbol{\mathcal{B}}_{m}(s)+\lambda\boldsymbol{\mathcal{B}}_{\mathcal{P}_{0}}(s)\\}$
$\displaystyle\quad+(\omega_{1i}-\omega_{0i})\\{\mathcal{A}_{m}(s)+\lambda\mathcal{P}_{1}(s)\mathcal{P}_{0}(s)\\}+\mathcal{P}_{0}(s)\mathcal{U}_{\lambda}(\textbf{D}_{i}).$

To derive the asymptotic distribution for $\widehat{\mathrm{PTE}},$ observe
that

$\displaystyle\widehat{\mathrm{PTE}}$ $\displaystyle-\mathrm{PTE}$
$\displaystyle=\int\left\\{\widehat{g}(s)-g_{\mathrm{opt}}(s)\right\\}d\left\\{\widehat{F}_{1}(s)-\widehat{F}_{0}(s)\right\\}+\int
g_{\mathrm{opt}}(s)d\left\\{\widehat{F}_{1}(s)-\widehat{F}_{0}(s)\right\\}-\mathrm{PTE}$
$\displaystyle=\frac{1}{\Delta}n^{-1}\sum_{i=1}^{n}\bigg{[}(\mathcal{G}_{m,i}(S_{i})+\lambda\mathcal{G}_{\mathcal{P}_{0},i}(S_{i}))(f_{1}(S_{i})-f_{0}(S_{i}))+(\omega_{1i}-\omega_{0i})g_{\mathrm{opt}}(S_{i})$
$\displaystyle\quad+\int
g_{\mathrm{opt}}(s)\\{\mathbf{f}^{\prime}_{1}(s)-\mathbf{f}^{\prime}_{0}(s)\\}\boldsymbol{\mathcal{U}}_{\alpha,i}ds-\mathrm{PTE}+\int\bigg{\\{}\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\\{\boldsymbol{\mathcal{B}}_{m}(s)+\lambda\boldsymbol{\mathcal{B}}_{\mathcal{P}_{0}}(s)\\}$
$\displaystyle\quad+(\omega_{1i}-\omega_{0i})\\{\mathcal{A}_{m}(s)+\lambda\mathcal{P}_{1}(s)\mathcal{P}_{0}(s)\\}+\mathcal{P}_{0}(s)\mathcal{U}_{\lambda}(\textbf{D}_{i})\bigg{\\}}(f_{1}(s)-f_{0}(s))ds\bigg{]}$
$\displaystyle\quad-\frac{\mathrm{PTE}}{\Delta}n^{-1}\sum_{i=1}^{n}[\omega_{1i}(Y_{i}-\mu_{1})-\omega_{0i}(Y_{i}-\mu_{0})]+o_{p}(n^{-1/2})$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\mathcal{U}_{\mathrm{PTE}}(\textbf{D}_{i})+o_{p}(n^{-1/2}),$

where

$\displaystyle\mathcal{U}_{\mathrm{PTE}}($ $\displaystyle\textbf{D}_{i})$
$\displaystyle=\frac{1}{\Delta}\bigg{[}(\mathcal{G}_{m,i}(S_{i})+\lambda\mathcal{G}_{\mathcal{P}_{0},i}(S_{i}))(f_{1}(S_{i})-f_{0}(S_{i}))+(\omega_{1i}-\omega_{0i})g_{\mathrm{opt}}(S_{i})$
$\displaystyle\quad+\int
g_{\mathrm{opt}}(s)\\{\mathbf{f}^{\prime}_{1}(s)-\mathbf{f}^{\prime}_{0}(s)\\}\boldsymbol{\mathcal{U}}_{\alpha,i}ds-\mathrm{PTE}+\int\bigg{\\{}\boldsymbol{\mathcal{U}}_{\alpha,i}^{\scriptscriptstyle\sf
T}\\{\boldsymbol{\mathcal{B}}_{m}(s)+\lambda\boldsymbol{\mathcal{B}}_{\mathcal{P}_{0}}(s)\\}$
$\displaystyle\quad+(\omega_{1i}-\omega_{0i})\\{\mathcal{A}_{m}(s)+\lambda\mathcal{P}_{1}(s)\mathcal{P}_{0}(s)\\}+\mathcal{P}_{0}(s)\mathcal{U}_{\lambda}(\textbf{D}_{i})\bigg{\\}}(f_{1}(s)-f_{0}(s))ds\bigg{]}$
$\displaystyle\quad-\frac{\mathrm{PTE}}{\Delta}[\omega_{1i}(Y_{i}-\mu_{1})-\omega_{0i}(Y_{i}-\mu_{0})].$

Therefore, by the central limit theorem,
$\sqrt{n}\left(\widehat{\mathrm{PTE}}-\mathrm{PTE}\right)$ converges in
distribution to a normal with mean zero and variance
$\sigma^{2}=E\left\\{\mathcal{U}_{\mathrm{PTE}}\left(\mathbf{D}_{i}\right)^{2}\right\\}$.

## Appendix 2

### Double Robustness

In this section, we prove that our proposed DR estimators are consistent when
either the PS model or the OR models are correctly specified. Recall that we
proposed the augmented IPW estimators for $m_{a}(s)$ and $f_{a}(s)$,

$\displaystyle\widehat{\mathcal{M}}_{\scriptscriptstyle a,\sf DR}(s)$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\left\\{K_{h}(S_{i}-s)Y_{i}\widehat{\omega}_{ai}-(\widehat{\omega}_{ai}-1)\widehat{\psi}^{\scriptscriptstyle{\dagger}}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{X}_{i})\widehat{\psi}^{\scriptscriptstyle{\dagger}}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{X}_{i})\right\\},$ (7a)
$\displaystyle\widehat{f}_{\scriptscriptstyle a,\sf DR}(s)$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\left\\{K_{h}(S_{i}-s)\widehat{\omega}_{ai}-(\widehat{\omega}_{ai}-1)\widehat{\psi}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{X}_{i})\right\\}$ (7b)
$\displaystyle\widehat{m}_{\scriptscriptstyle a,\sf DR}(s)$
$\displaystyle=\frac{\widehat{\mathcal{M}}_{\scriptscriptstyle a,\sf
DR}(s)}{\widehat{f}_{\scriptscriptstyle a,\sf DR}(s)}$ (7c)

respectively, where $h=O(n^{-\nu})$ with $\nu\in(1/4,1/2)$,
$\widehat{\psi}_{a,\scriptscriptstyle\sf m}(s;\mathbf{x})$ and
$\widehat{\psi}_{a,\scriptscriptstyle\sf f}(s;\mathbf{x})$ are estimators for
the conditional mean $\psi_{a,\scriptscriptstyle\sf
m}^{\scriptscriptstyle{\dagger}}(s;\mathbf{x})$ and the conditional density
$\psi_{a,\scriptscriptstyle\sf
f}^{\scriptscriptstyle{\dagger}}(s;\mathbf{x})$, respectively.

We now show that the estimators are consistent if either
$\sup_{\mathbf{x}}|\pi_{a}(\mathbf{x};\widehat{\boldsymbol{\alpha}})-\pi_{a}(\mathbf{x})|\to
0$ in probability or
$\sup_{\mathbf{x},s}\\{|\widehat{\psi}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{x})-\psi_{a,\scriptscriptstyle\sf
m}^{\scriptscriptstyle{\dagger}}(s;\mathbf{x})|+|\widehat{\psi}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{x})-\psi_{a,\scriptscriptstyle\sf
f}^{\scriptscriptstyle{\dagger}}(s;\mathbf{x})|\\}\to 0$ in probability. Let
$\bar{\boldsymbol{\alpha}}$, $\bar{\psi}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{x})$, $\bar{\psi}_{a,\scriptscriptstyle\sf f}(s;\mathbf{x})$
denote the respective limits of $\widehat{\boldsymbol{\alpha}}$,
$\widehat{\psi}_{a,\scriptscriptstyle\sf m}(s;\mathbf{x})$ and
$\widehat{\psi}_{a,\scriptscriptstyle\sf f}(s;\mathbf{x})$ under possible mis-
specification of their respective models,
$\bar{\pi}_{a}(\mathbf{x})=\pi_{a}(\mathbf{x};\bar{\boldsymbol{\alpha}})$, and
$\bar{\omega}_{ai}=I(A_{i}=a)/\bar{\pi}_{a}(\mathbf{X}_{i})$. Regardless of
the adequacy of the models, by the central limit theorem and convergence of
kernel smoothed estimators (Pagan & Ullah, 1999), we have that
$\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}}=O_{p}(n^{-\frac{1}{2}})$
and $\sup_{s,\mathbf{x}}|\widehat{\psi}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{x})-\bar{\psi}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{x})|+\sup_{s,\mathbf{x}}|\widehat{\psi}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{x})-\bar{\psi}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{x})|=o_{p}(1)$.

When the PS model is correctly specified,
$\sup_{s}|\widehat{m}_{a}(s)-m_{a}(s)|+\sup_{s}|\widehat{f}_{a}(s)-f_{a}(s)|\to
0$ in probability as shown in Appendix 1. In addition, the augmentation terms

$n^{-1}\sum_{i=1}^{n}(\widehat{\omega}_{ai}-1)\widehat{\psi}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{X}_{i})=n^{-1}\sum_{i=1}^{n}(\omega_{ai}-1)\bar{\psi}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{X}_{i})+O_{p}(\|\widehat{\boldsymbol{\alpha}}-\boldsymbol{\alpha}_{0}\|_{2})$

and

$n^{-1}\sum_{i=1}^{n}(\widehat{\omega}_{ai}-1)\widehat{\psi}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{X}_{i})=n^{-1}\sum_{i=1}^{n}(\omega_{ai}-1)\bar{\psi}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{X}_{i})+O_{p}(\|\widehat{\boldsymbol{\alpha}}-\boldsymbol{\alpha}_{0}\|_{2})$

also converge to 0 in probability, regardless of the adequacy of the OR
models. Therefore, under the correct specification of the PS model,
$\sup_{s}|\widehat{m}_{\scriptscriptstyle a,\sf
DR}(s)-m_{a}(s)|+\sup_{s}|\widehat{f}_{\scriptscriptstyle a,\sf
DR}(s)-f_{a}(s)|\to 0$ in probability.

We next establish the consistency of the DR estimators when the PS model may
be mis-specified but the OR models are correctly specified. First consider
$\widehat{f}_{\scriptscriptstyle a,\sf DR}(s)$, which can be written as

$\displaystyle\widehat{f}_{\scriptscriptstyle a,\sf DR}(s)$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\bigg{[}\\{K_{h}(S_{i}-s)-\psi_{a,\scriptscriptstyle\sf
f}^{\scriptscriptstyle{\dagger}}(s;\mathbf{X}_{i})\\}\widehat{\omega}_{ai}$
$\displaystyle\quad-(\widehat{\omega}_{ai}-1)\\{\widehat{\psi}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{X}_{i})-\psi_{a,\scriptscriptstyle\sf
f}^{\scriptscriptstyle{\dagger}}(s;\mathbf{X}_{i})\\}+\psi_{a,\scriptscriptstyle\sf
f}^{\scriptscriptstyle{\dagger}}(s;\mathbf{X}_{i})\bigg{]}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\bigg{[}\\{K_{h}(S_{i}-s)-\psi_{a,\scriptscriptstyle\sf
f}^{\scriptscriptstyle{\dagger}}(s;\mathbf{X}_{i})\\}\widehat{\omega}_{ai}$
$\displaystyle\quad-(\widehat{\omega}_{ai}-1)\\{\widehat{\psi}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{X}_{i})-\psi_{a,\scriptscriptstyle\sf
f}^{\scriptscriptstyle{\dagger}}(s;\mathbf{X}_{i})\\}\bigg{]}+f_{a}(s)+O_{p}(n^{-\frac{1}{2}}),$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\varepsilon_{\scriptscriptstyle a,\sf
DR}(s;\mathbf{D}_{i})+f_{a}(s)-n^{-1}\sum_{i=1}^{n}(\bar{\omega}_{ai}-1)\\{\widehat{\psi}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{X}_{i})-\psi_{a,\scriptscriptstyle\sf
f}^{\scriptscriptstyle{\dagger}}(s;\mathbf{X}_{i})\\}$
$\displaystyle\quad+O_{p}(n^{-\frac{1}{2}}),$

where $\varepsilon_{a,\scriptscriptstyle\sf
f}(s;\mathbf{D}_{i})=\\{K_{h}(S_{i}-s)-\psi_{a,\scriptscriptstyle\sf
f}^{\scriptscriptstyle{\dagger}}(s;\mathbf{X}_{i})\\}\bar{\omega}_{ai}$. It
follows from the uniform convergence of kernel smoothed estimators (Pagan &
Ullah, 1999) that

$\sup_{s}|n^{-1}\sum_{i=1}^{n}\varepsilon_{\scriptscriptstyle a,\sf
DR}(s;\mathbf{D}_{i})-E\\{\varepsilon_{\scriptscriptstyle a,\sf
DR}(s;\mathbf{D}_{i})\\}|=o_{p}(1)$

and

$E\\{\varepsilon_{\scriptscriptstyle a,\sf
DR}(s;\mathbf{D}_{i})\\}=E\left[\bar{\omega}_{ai}\left\\{\int
K_{h}(S-s)\psi_{a,\scriptscriptstyle\sf
f}^{\scriptscriptstyle{\dagger}}(S;\mathbf{X}_{i})dS-\psi_{a,\scriptscriptstyle\sf
f}^{\scriptscriptstyle{\dagger}}(s;\mathbf{X}_{i})\right\\}\right]=O(h^{2}).$

This together with
$\sup_{s,\mathbf{x}}|\widehat{\psi}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{x})-\psi^{\scriptscriptstyle{\dagger}}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{x})|=o_{p}(1)$ implies that
$\sup_{s}|\widehat{f}_{\scriptscriptstyle a,\sf DR}(s)-f_{a}(s)|=o_{p}(1)$. We
have a similar consistency result for
$\widehat{\mathcal{M}}_{\scriptscriptstyle a,\sf DR}(s)$, where

$\displaystyle\widehat{\mathcal{M}}_{\scriptscriptstyle a,\sf DR}(s)$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\left\\{K_{h}(S_{i}-s)Y_{i}\bar{\omega}_{ai}-(\bar{\omega}_{ai}-1)\psi^{\scriptscriptstyle{\dagger}}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{X}_{i})\psi^{\scriptscriptstyle{\dagger}}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{X}_{i})\right\\}+o_{p}(1)$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\bigg{[}\varepsilon_{a,\scriptscriptstyle\sf
m}(s;\mathbf{D}_{i})+\left\\{K_{h}(S_{i}-s)\psi_{a,\scriptscriptstyle\sf
m}^{\scriptscriptstyle{\dagger}}(s;\mathbf{X}_{i})-\psi^{\scriptscriptstyle{\dagger}}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{X}_{i})\psi^{\scriptscriptstyle{\dagger}}_{a,\scriptscriptstyle\sf
f}(s;\mathbf{X}_{i})\right\\}\bar{\omega}_{ai}\bigg{]}$
$\displaystyle\quad+m_{a}(s)+o_{p}(1),$

and $\varepsilon_{a,\scriptscriptstyle\sf
m}(s;\mathbf{D}_{i})=K_{h}(S_{i}-s)\\{Y_{i}-\psi_{a,\scriptscriptstyle\sf
m}^{\scriptscriptstyle{\dagger}}(s;\mathbf{X}_{i})\\}\bar{\omega}_{ai}$.
Following the convergence of $\widehat{f}_{\scriptscriptstyle a,\sf DR}(s)\to
f_{a}(s)$, $\widehat{\mathcal{M}}_{\scriptscriptstyle a,\sf DR}(s)\to
m_{a}(s)f_{a}(s)$, $\widehat{\psi}_{a,\scriptscriptstyle\sf
m}(s;\mathbf{x})\to\psi_{a,\scriptscriptstyle\sf
m}^{\scriptscriptstyle{\dagger}}(s;\mathbf{x})$ and
$\widehat{\omega}_{ai}\to\bar{\omega}_{ai}$, we arrive at the consistency of
$\widehat{m}_{\scriptscriptstyle a,\sf DR}(s)$ to $m_{a}(s)$ when the PS model
may be mis-specified but the OR models are correctly specified.

Thus, we get the double robustness properties for
$\widehat{f}_{\scriptscriptstyle a,\sf DR}(s)$ and
$\widehat{m}_{\scriptscriptstyle a,\sf DR}(s)$.

Since all remaining estimators relevant to $\widehat{g}_{\scriptscriptstyle\sf
DR}(s)$ are plug-in estimators that are derived based on
$\widehat{m}_{\scriptscriptstyle a,\sf DR}(s)$ and
$\widehat{f}_{\scriptscriptstyle a,\sf DR}(s)$, we can conclude the double
robustness of $\widehat{g}_{\scriptscriptstyle\sf DR}(s)$ for
$g_{\scriptscriptstyle\sf opt}(s)$.

Finally, the PTE will be doubly robust by standard arguments for the
conditional mean estimators (Robins et al., 1994), where we construct a plug-
in estimator for $\Delta_{g_{\scriptscriptstyle\sf opt}}$ as
$\widehat{\Delta}_{\widehat{g},\scriptscriptstyle\sf
DR}=\widehat{\mu}_{1,\widehat{g},\scriptscriptstyle\sf
DR}-\widehat{\mu}_{0,\widehat{g},\scriptscriptstyle\sf DR},$ where

$\widehat{\mu}_{a,g,DR}=n_{a}^{-1}\sum_{i:A_{i}=a}\left\\{\frac{g(S_{i})}{\hat{\pi}_{a}(\mathbf{X}_{i})}-\frac{I(A_{i}=a)-\hat{\pi}_{a}(\mathbf{X}_{i})}{\hat{\pi}_{a}(\mathbf{X}_{i})}\hat{\zeta}_{a,g}(\mathbf{X}_{i})\right\\},$

where $\hat{\zeta}_{a,g}(\mathbf{x})$ is an estimator for
$\zeta_{a,g}(\mathbf{x})=E(g(S_{i}^{(a)})\mid\mathbf{X}_{i}=\mathbf{x})=E(g(S_{i})\mid
A_{i}=a,\mathbf{X}_{i}=\mathbf{x})$, and $n_{a}=\sum_{i=1}^{n}I(A_{i}=a)$,
$a=0,1$. Similarly, we define $\widehat{\Delta}_{\scriptscriptstyle\sf
DR}=\widehat{\mu}_{\scriptscriptstyle 1,\sf
DR}-\widehat{\mu}_{\scriptscriptstyle 0,\sf DR}$, where

$\widehat{\mu}_{\scriptscriptstyle a,\sf
DR}=n_{a}^{-1}\sum_{i:A_{i}=a}\left\\{\frac{Y_{i}}{\widehat{\pi}_{a}(\mathbf{X}_{i})}-\frac{I(A_{i}=a)-\widehat{\pi}_{a}(\mathbf{X}_{i})}{\widehat{\pi}_{a}(\mathbf{X}_{i})}\hat{\zeta}_{a}(\mathbf{X}_{i})\right\\},$

where $\hat{\zeta}_{a}(\mathbf{x})$ is an estimator for
$\zeta_{a}(\mathbf{x})=E(Y_{i}^{(a)}\mid\mathbf{X}_{i}=\mathbf{x})=E(Y_{i}\mid
A_{i}=a,\mathbf{X}_{i}=\mathbf{x})$.