11affiliationtext: Department of Epidemiology, Harvard T.H. Chan School of
Public Health, Boston, MA. 22affiliationtext: CAUSALab, Harvard T.H. Chan
School of Public Health, Boston, MA. 33affiliationtext: Department of
Biostatistics, Harvard T.H. Chan School of Public Health, Boston, MA.
44affiliationtext: Department of Biostatistics, Brown University School of
Public Health, Providence, RI.

# Assessing model performance for counterfactual predictions

Christopher B. Boyer email<EMAIL_ADDRESS>Issa J. Dahabreh  Jon A.
Steingrimsson

Counterfactual prediction methods are required when a model will be deployed
in a setting where treatment policies differ from the setting where the model
was developed, or when the prediction question is explicitly counterfactual.
However, estimating and evaluating counterfactual prediction models is
challenging because one does not observe the full set of potential outcomes
for all individuals. Here, we discuss how to tailor a model to a
counterfactual estimand, how to assess the model’s performance, and how to
perform model and tuning parameter selection. We also provide identifiability
results for measures of performance for a potentially misspecified
counterfactual prediction model based on training and test data from the same
(factual) source population. Last, we illustrate the methods using simulation
and apply them to the task of developing a statin-naïve risk prediction model
for cardiovascular disease.

Keywords: causal inference, prediction model, treatment drop-in,
transportability, model performance, machine learning

##  1 Introduction

Many common tasks in prediction modeling can be viewed as a form of
counterfactual prediction. First, prediction models are often deployed in
settings that differ from those in which they are trained and one of the ways
settings may differ is in terms of treatment patterns after baseline [1]. For
example, a model fit in a setting where 5% are treated over the follow up
period may not produce valid predictions in a setting where 50% are treated
and vice versa. Even when models are deployed in the same setting, treatment
policies may change over time, affecting who is likely to be treated and
leading to problems of “domain adaption” or “dataset shift” [2, 3]. These
differences between training and deployment settings can cause the performance
of models to degrade, particularly when, as is often the case, model
predictors are themselves correlated with, or direct determinants of,
treatment [4, 5, 6].

When faced with a change in treatment policy, ideally, one would re-train the
model using data reflecting the current patterns of treatment in the target
population where the model will be applied. Collecting the necessary data to
re-train the model, however, may be inordinately expensive or time consuming.
Absent sufficient resources or as a stop gap, one might consider using
existing training data to tailor the model to target the expected outcome that
would be observed were treatment to be administered as in the deployment
setting [7]. Alternatively, one might want to estimate how poorly the existing
model is likely to perform in the deployment setting, to determine whether
additional data collection efforts are worthwhile. In either case, the
research questions pertain to counterfactual predictions.

Second, there are instances in which the prediction question itself is
inherently counterfactual. For instance, clinicians may use a model to counsel
patients about their risk of disease if they were to remain untreated or to
compare risks under alternative treatment strategies [8, 9]. In some
circumstances, such as when data from a randomized trial are used to model
conditional average treatment effects, this can be done without additional
formalism [10, 11, 12, 13, 14]. In contrast, when training data are obtained
from an observational study where treatment initiation over follow up is not
strictly controlled by the investigators, the predictions most relevant to
decision-making are counterfactual [15, 9].

In both instances, when treatment policies change across settings or when the
prediction question is inherently counterfactual, we need methods for
tailoring models to answer counterfactual prediction questions, when data on
the full set of potential outcomes is not available. We also need performance
measures that evaluate model performance in these new environments without
assuming the prediction model is correctly specified. Here, we examine the
conditions under which tailoring a model for counterfactual prediction is
possible using training data alone. Under similar conditions, we also show
that the counterfactual performance of the model may be estimated
independently from the method used to fit the model and may be evaluated even
if the model is misspecified. Performance measures can therefore be used to
differentiate between better and worse-performing models or to quantify how
badly a model is likely to perform in a new setting.

##  2 Set up and notation

Let $Y$ be the outcome of interest, $X$ a vector of baseline covariates, and
$A$ an indicator of treatment over the follow up period. We assume that data
are obtained from a simple random sample from a target population
$\\{(X_{i},A_{i},Y_{i})\\}_{i=1}^{n}$ in which treatment is initiated
according to pattern $f^{obs}(A|X)$. Covariates in $X$ include a set
sufficient to control confounding of the treatment-outcome relationship ($L$)
as well as additional predictors of the outcome ($P$). Our goal is to build a
prediction model for $Y$ using only covariates $X^{*}$ which are a subset of
$X$, chosen based on their availability and predictive potential rather than
on whether they can help control confounding (note that $X^{*}$ can include
components of both $L$ and $P$). An example causal directed acyclic graph for
this process is shown in Figure 1.

$L$$A$$Y$$P$$U$ Figure 1: Example causal directed acyclic graph for prediction
in a setting with a single time fixed treatment $A$ over follow up.

The data are randomly split into a training set of $n_{train}$ observations
and a test set of $n_{test}$ observations, with $n=n_{train}+n_{test}$. Let
$D_{train}$ and $D_{test}$ be indicators of whether an observation is in the
training set or test set respectively. As is customary, we use the training
set to build a prediction model for the expected outcome conditional on
covariates $\operatorname{E}[Y|X^{*}]$ and, then use the test set to evaluate
model performance. Our results, however, apply equally to existing prediction
models or (directly measured) biomarkers, in which case splitting is not
necessary. Let $\mu_{\beta}(X^{*})$ be a parametric model, indexed by
parameter $\beta$, and $\mu_{\widehat{\beta}}(X^{*})$ be the “fitted” model
using parameter estimates $\widehat{\beta}$. We allow for the possibility that
model $\mu_{\beta}(X^{*})$ is misspecified. For a particular estimand such as
$\operatorname{E}[Y|X^{*}]$, a model is correctly specified if there exists
$\beta_{0}\in\mathcal{B}$, where $\mathcal{B}$ is the parameter space of
$\beta$, such that $\mu_{\beta_{0}}(X^{*})=\operatorname{E}[Y|X^{*}]$ and the
model is misspecified if no such $\beta_{0}$ exists. We use $f(\cdot)$
generically to denote densities.

To define counterfactual estimands of interest, let $Y^{a}$ be the potential
outcome under an intervention which sets treatment $A$ to $a$. To keep our
notation simple, here we focus on the case that treatment is time-fixed – it
is initiated after baseline and sustained throughout or its effect is
independent of duration. In Section B of the Appendix, we extend our setup to
handle time-varying treatments. We also assume interventions on $A$ are
deterministic, but extend to random and dynamic interventions in Section A.4
of the Appendix. These may be more relevant for some real world prediction
tasks, for instance, when evaluating the performance of models in settings
where treatment differences are stochastic. Finally, we also assume that there
is no loss to follow up, but note that our results can be extended by
recognizing censoring can be accommodated using the methods we proposed for
time-varying treatment [16] for a hypothetical intervention to prevent loss-
to-followup.

##  3 Estimands for assessing model performance

Our goal is to train a model for counterfactual prediction and determine the
model’s performance. We posit a parametric prediction model
$\mu_{\beta}(X^{*})$ for the expected potential outcome under intervention
$A=a$ conditional on covariates $\operatorname{E}[Y^{a}|X^{*}]$, which we wish
to estimate from the training dataset. The model may be misspecified, such as
if, for instance, it is a model trained on another target such as the expected
(factual) outcome in the source population $\operatorname{E}[Y|X^{*}]$.

To evaluate prediction model performance, a number of measures have been
proposed in the prediction literature [17, 18, 19]. These measures generally
compare fitted predictions $\mu_{\widehat{\beta}}(X^{*})$ and the observed
outcomes $Y$ via a distance function and associated sample statistic. For
counterfactual predictions, however, model evaluation is more challenging
because the potential outcome $Y^{a}$ is not directly observed. Yet, as we
will show, under certain conditions the expected value of the performance
measure may still be identified from the observed data in the test set. An
example performance measure of interest is

$\psi(a)=\operatorname{E}[(Y^{a}-\mu_{\widehat{\beta}}(X^{*}))^{2}\mid
D_{test}=1]$

where the squared error loss $(Y^{a}-\mu_{\widehat{\beta}}(X^{*}))^{2}$
quantifies the discrepancy between the potential outcome under treatment level
$A=a$ and the model prediction $\mu_{\widehat{\beta}}(X^{*})$ in terms of the
squared difference. In the main text, we focus on the mean squared error as
the performance measure of interest $\psi(a)$. In the Appendix, we extend our
results to the case where the measure is any member of a generic class of
counterfactual loss functions $L(Y^{a},\mu_{\widehat{\beta}}(X^{*}))$ (Section
A.2) as well as more complex risk-based measures such as area under the
receiver operating characteristics curve (AUC, Section E), which depends on
paired observations, and the calibration curve (Section D), which is a
functional. Importantly, all the measures considered are identifiable without
assuming that the model $\mu_{\widehat{\beta}}(X^{*})$ is correctly specified.

##  4 Identifiability conditions

We will assume the following identifiability assumptions which have been
described in more detail elsewhere [16, 20, 21].

  1. A1.

Consistency. If $A=a$, then $Y^{a}=Y$

  2. A2.

Conditional exchangeability. $Y^{a}\perp\\!\\!\\!\perp A\mid X$

  3. A3.

Positivity. For all $x$ with positive density, i.e. $f_{X}(x)>0$, $\Pr[A=a\mid
X=x]>0$

Consistency implies that observed outcomes among those with $A=a$ reflect
potential outcomes under corresponding level of treatment. It would be
violated if, for instance, there were multiple “hidden” versions of the
treatment under consideration [22]. The exchangeability condition stipulates
that treatment is conditionally independent of the potential outcome given
covariates $X$. Finally, the positivity condition implies that there is a
positive probability of observed treatment level $A=a$ in all strata of $X$.

##  5 Tailoring a model for counterfactual predictions

As we show in Appendix Section A.1, under the conditions above
$\operatorname{E}[Y^{a}|X^{*}]$ is identified by

$\phi_{a}(X^{*})\equiv\operatorname{E}[\operatorname{E}[Y\mid
X,A=a,D_{train}=1]\mid X^{*},D_{train}=1]$ (1)

or, equivalently, using an inverse probability weighted expression

$\phi_{a}(X^{*})=\operatorname{E}\left[\frac{I(A=a)}{\Pr[A=a\mid
X,D_{train}=1]}Y\Big{|}X^{*},D_{train}=1\right]$ (2)

The two expressions for $\phi_{a}(X^{*})$ suggest possible approaches for
tailoring the model for counterfactual predictions using only the training
data.

One approach, based on equation 1, is to subset to participants with
corresponding treatment level $A=a$ in the training data and fit a model
$\mu_{\beta}(X)$ for the observed $Y$ conditional $X$, i.e.
$\operatorname{E}[Y|X,A=a]=\mu_{\beta}(X)$. Then when the desired predictors
$X^{*}$ are a subset of $X$, the covariates sufficient to ensure
exchangeability, predictions are marginalized (standardized) over the
covariates in $X$ that are not in $X^{*}$. When the dimension of $X^{*}$ is
small, this can be done nonparametrically; however, when $X^{*}$ is higher-
dimensional, an additional modeling step will be required either (1) modeling
the estimated $\widehat{\mu}_{\widehat{\beta}}(X)$ as a function of $X^{*}$,
i.e. $\operatorname{E}[\widehat{\mu}_{\widehat{\beta}}(X)|X^{*}]$ or (2)
modeling the conditional density of $X$ given $X^{*}$, i.e. $f(X|X^{*})$. The
resulting predictions will be consistent for $\operatorname{E}[Y^{a}|X^{*}]$
provided all models are correctly specified. An alternative suggested in [7],
would be to simulate samples from the model $\mu_{\beta}(X)$ and fit a second
stage model $\mu^{*}_{\beta}(X^{*})$ using only the subset of the predictors
of interest. This is similar to the first standardization method above and has
the advantage that the second stage model may be developed using an outcome
with same support as the original $Y$.

A second approach, based on 2, is to fit a weighted model
$\mu_{\beta}(X^{*})$, using for instance weighted maximum likelihood, with
weights equal to the probability of receiving treatment level $A=a$
conditional on covariates $X$ necessary to ensure exchangeability, i.e.,
sample analogs of $W(a)=\frac{I(A=a)}{\Pr[A=a\mid X,D_{train}=1]}$. This is
the basis for previously proposed methods for counterfactual prediction based
on inverse probability of treatment weighting [23, 1]. Note that, as before,
it is possible to specify a subset of predictors $X^{*}$ used in the
prediction model $\mu_{\beta}(X^{*})$ as compared to the full set of
covariates $X$ required for exchangeability which are only necessary for
defining the weights $W$, however a second marginalization or simulation step
is not required. This means tailoring the model for counterfactual predictions
using the weighting approach can be accomplished more easily using off-the-
shelf software.

##  6 Assessing model performance

Using the same identifiability conditions, in Appendix Section A.2, we show
the model performance measure $\psi(a)$ is identifiable using data from the
test set through the expression

$\psi(a)\equiv\operatorname{E}\left[\operatorname{E}[(Y-\mu_{\widehat{\beta}}(X^{*}))^{2}\mid
X,A=a,D_{test}=1]\mid D_{test}=1\right]$ (3)

or, equivalently, using an inverse probability weighted expression,

$\psi(a)=\operatorname{E}\left[\frac{I(A=a)}{\Pr[A=a\mid
X,D_{test}=1]}(Y-\mu_{\widehat{\beta}}(X^{*}))^{2}\mid D_{test}=1\right]$ (4)

regardless of whether the model $\mu_{\widehat{\beta}}(X^{*})$ has been
tailored to target $\operatorname{E}[Y^{a}\mid X]$ or is correctly specified
in general.

As before, the two expression suggest two different approaches for the
estimation of model performance using the test data alone.

First, using the sample analog of expression (3), an estimator of the target
MSE is

$\widehat{\psi}_{CL}(a)=\frac{1}{n_{test}}\sum_{i=1}^{n}I(D_{test,i}=1)\widehat{h}_{a}(X_{i})$
(5)

where $\widehat{h}_{a}(X)$ is an estimator for the conditional loss
$\operatorname{E}[(Y-\mu_{\widehat{\beta}}(X^{*}))^{2}\mid X,A=a,D_{test}=1]$
and from the perspective of estimating $\psi(a)$ may be considered a nuisance
function. To keep notation simple, we suppress the dependency of
$\widehat{h}_{a}(X)$ on $\mu_{\widehat{\beta}}$. When the dimension of $X$ is
small it may be possible to use sample analogs of $\widehat{h}_{a}(X)$. In
almost all practical cases though some form of modeling will be required; in
these cases, $\widehat{\psi}_{CL}(a)$ is a consistent estimator for $\psi(a)$
as long as the model for $\widehat{h}_{a}(X)$ is correctly specified.

Next, using the sample analog of expression (4), an alternative weight-based
estimator of the target MSE is

$\widehat{\psi}_{IPW}(a)=\frac{1}{n_{test}}\sum_{i=1}^{n}\frac{I(A_{i}=a,D_{test,i}=1)}{\widehat{e}_{a}(X_{i})}(Y_{i}-\mu_{\widehat{\beta}}(X^{*}_{i}))^{2}$
(6)

where $\widehat{e}_{a}(X)$ is another nuisance function estimating the
probability of receiving treatment level $A=a$ conditional on $X$, i.e.
$\Pr[A=a\mid X,D_{test}=1]$. Again, when the dimension of $X$ is small it may
be possible to use the sample analog of $\widehat{e}_{a}(X)$, but in most
cases, it will have to be modeled. In these cases, $\widehat{\psi}_{IPW}(a)$
is a consistent estimator of $\psi(a)$ as long as the model for
$\widehat{e}_{a}(X)$ is correctly specified.

The conditional loss estimator (5) relies on correctly specifying the model
for the conditional loss and the weighting estimator (6) relies on correctly
specifying the model for the probability of treatment. In some settings, one
estimator may be preferred over the other: when more is known about the
mechanism for “assigning” treatment or when the outcome is very rare, the
weighting estimator may be preferred [24, 25]; when the process that gives
rise to the outcomes is well understood the conditional loss estimator may be
preferred. In practice, however, both models may be difficult to specify
correctly. Using data-adaptive and more flexible machine learning estimators
for estimation of these nuisance models offers the possibility of capturing
more complex data generation processes. These data-adaptive estimators
generally have slower rates of convergence than the $\sqrt{n}$ rates of
parametric models and therefore will not yield asymptotically valid confidence
intervals [26]. To address this challenge, we can use a doubly-robust
estimator which combines models for the conditional loss $\widehat{h}_{a}(X)$
and the probability of treatment $\widehat{e}_{a}(X)$, such as

$\widehat{\psi}_{DR}(a)=\frac{1}{n_{test}}\sum_{i=1}^{n}I(D_{test,i}=1)\left[\widehat{h}_{a}(X_{i})+\frac{I(A_{i}=a)}{\widehat{e}_{a}(X_{i})}\left\\{(Y_{i}-\mu_{\widehat{\beta}}(X^{*}_{i}))^{2}-\widehat{h}_{a}(X_{i})\right\\}\right]$
(7)

As we show in Appendix section C, under mild regularity conditions [27], this
estimator will be consistent if one of $\widehat{h}_{a}(X)$ and
$\widehat{e}_{a}(X)$ is correctly specified. It also permits the use of
machine learning or data-adaptive estimators that converge at rate slower than
$\sqrt{n}$, thus allowing for more flexible estimation of the nuisance
functions. This is due to the fact that the empirical process terms governing
the convergence of $\widehat{\psi}_{DR}$ involve a product of the estimation
errors for $\widehat{h}_{a}(X)$ and $\widehat{e}_{a}(X)$ which converge under
the weaker condition that only the combined rate of convergence for both
nuisance functions is at least $\sqrt{n}$ [26].

##  7 Model and tuning parameter selection

Up to this point, we have assumed that $\mu_{\beta}(X^{*})$ is a pre-specified
parametric model and ignored any form of model selection (e.g. variable or
other specification search) or data-adaptive tuning parameter selection, which
may be the case when using an existing (validated) model. In many cases,
however, analysts have to select between multiple models or perform a data-
adaptive search through a parameter space for tuning parameter selection when
developing a prediction model [19]. To avoid overfitting, analysts typically
use methods such as cross-validation or the bootstrap to perform selection.
These techniques rely on optimizing some measure of model performance, such as
the MSE.

When performing model or tuning parameter selection for counterfactual
prediction, the results from the previous sections suggest that the model
performance measure should be targeted to the counterfactual performance if
the intervention of interest were universally applied. For example, when using
cross-validation for model selection the analyst splits the data into $K$
mutually exclusive “folds” and fits the candidate models using $K-1$ of the
folds and estimates the performance of each in the held out fold. This process
is repeated $K$ times where each fold is left out once. The final performance
estimate is the average of the $K$ estimates and the model with best overall
performance is selected (or, alternatively, the tuning parameter with the best
performance). When performing counterfactual prediction, at each stage in the
procedure the analyst should use modified performance measures such as those
in section 6 above. Failure to do so, can lead to sub-optimal selection with
respect to the counterfactual prediction of interest.

##  8 Simulation experiments

We performed two Monte Carlo simulation experiments to illustrate (i) the
benefits of tailoring models to the correct counterfactual estimand of
interest, (ii) the potential for bias when using naïve estimators of model
performance such as the MSE, (iii) the importance of correct specification of
the nuisance models when estimating counterfactual performance, and (iv) the
properties of the doubly-robust estimator under misspecification of the
nuisance models. We adapt data generation processes previously used for
transporting models between settings under covariate shift [28, 29].

###  8.1 Experiment 1

We simulated treatment initiation at baseline based on the logistic model
$\Pr[A=1\mid X]=\operatorname{expit}(1.5-0.3X)$, where predictors $X$ are
drawn from $X\sim$ Uniform $(0,10)$. Under this model, about 50% initiate
treatment but those with higher values of $X$ are less likely to start
treatment than those with lower values of $X$. We then simulated the outcome
using the linear model $Y=1+X+0.5X^{2}-3A+\varepsilon$, where
$\varepsilon\sim\mathcal{N}(0,X)$. We set the total sample size to 1000 and
the data were randomly split in a 1:1 ratio into a training and a test set. In
this simulation, $X=X^{*}$. The full process may be written:

$\displaystyle X$ $\displaystyle\sim\text{Unif}(0,10)$ $\displaystyle A$
$\displaystyle\sim\text{Bernoulli}\\{\operatorname{expit}(-1.5+0.3X)\\}$
$\displaystyle Y$ $\displaystyle\sim\text{Normal}(1+X+0.5X^{2}-3A,X)$

Our goal was to estimate a model in the same population where, contrary to
fact, treatment was universally withheld, i.e. we targeted
$\operatorname{E}[Y^{a=0}\mid X]$. Note that, for simplicity, in this case
predictor variables are the same set necessary to control confounding
$X=X^{*}$. Under this data generating mechanism, the MSE under no treatment is
larger than the MSE under the natural course and identifiability conditions
1-3 are satisfied. We considered two specifications of prediction models
$\mu_{\beta}(X)$:

  1. 1.

a correctly specified linear regression model that included the main effects
of $X$ and $X^{2}$, i.e. $\mu_{\beta}(X)=\beta_{0}+\beta_{1}X+\beta_{2}X^{2}$.

  2. 2.

a misspecified linear regression model that only included the main effect of
$X$, i.e. $\mu_{\beta}(X)=\beta_{0}+\beta_{1}X$.

For each specification, we also considered two estimation strategies: one
using ordinary least squares regression (OLS) and ignoring treatment
initiation and the other using weighted least squares regression (WLS) where
the weights were equal to the inverse of the probability of being untreated.
As discussed above the latter specifically targets the counterfactual estimand
under no treatment. Finally, we considered two approaches for estimating the
performance of the models in the test set: a naïve estimate of the MSE using
observed outcome values, i.e.

$\widehat{\psi}_{Na\ddot{i}ve}(a)=\frac{1}{n_{test}}\sum_{i=1}^{n}I(D_{test,i}=1)(Y_{i}-\mu_{\widehat{\beta}}(X_{i}))^{2},$

and the inverse-probability weighted estimator $\widehat{\psi}_{IPW}(a)$ from
section 6. For the latter, we fit a correctly specified logistic regression
model for $e_{a}(X)$, i.e.
$e_{a}(X)=\operatorname{expit}(\alpha_{0}+\alpha_{1}X)$, in the test set to
estimate the weights. Lastly, we also calculated the “true” MSE under
intervention to withold treatment by generating test data under same process
as above but setting the treatment $A$ to $a=0$ for everyone and then
averaging across simulations.

Table 1: Simulation results for experiment 1 comparing the performance of OLS
and WLS models using the naïve and inverse probability weighting (IPW)
estimators of the MSE.

Model $\mu_{\beta}(X)$ |  Naïve |  IPW |  Truth
---|---|---|---
Correct |  |  |
OLS |  2.9 |  3.6 |  3.6
WLS |  5.5 |  1.0 |  1.0
Misspecified |  |  |
OLS |  16.8 |  17.5 |  17.5
WLS |  19.5 |  15.0 |  15.0

  * •

Correct and misspecified refers to the specification of the prediction model
$\mu_{\beta}(X)$. OLS = model estimation using ordinary least squares
regression (unweighted); WLS = model estimation using weighted least squares
regression with weights equal to the inverse probability of being untreated.
Results were averaged over 10,000 simulations. The true counterfactual MSE was
obtained using numerical methods.

Table 1 shows the results of the experiment based on 10,000 simulations. In
general, correctly specified models yielded smaller average MSE than
misspecified models. Comparing the performance of OLS and WLS estimation, when
using $\widehat{\psi}_{Na\ddot{i}ve}(a)$, the naïve estimator of the MSE, OLS
seemed to produce better predictions than WLS when correctly specified
(average MSE of 2.9 vs. 5.5) as well as when misspecified (average MSE of 16.8
vs. 19.5). In contrast, when using $\widehat{\psi}_{IPW}(a)$, the inverse-
probability weighted estimate of the MSE, WLS performed better than OLS both
when the model was correctly specified (average MSE of 1.0 vs. 3.6) and when
misspecified (average MSE of 15.0 vs. 17.5). For reference, in the last column
we show the true counterfactual MSE that would be obtained if one had access
to the potential outcomes (obtained via numerical methods). We found that the
average of the inverse probability weighted estimator across the simulations
was equivalent to this quantity for all specifications and for both OLS and
WLS estimation. This suggests that only the modified estimators of model
performance in section 6 are able to accurately estimate the counterfactual
performance of the model. Indeed, under this data generation process, if one
were to use the naïve estimator one might erroneously conclude that the OLS
model is the better choice.

###  8.2 Experiment 2

In the previous experiment, we assumed the nuisance models for the MSE were
correctly specified. We now consider estimation of model performance in the
more likely case that nuisance models are misspecified. Using the results from
Section E in the Appendix, we also estimate the area under the receiver
operating characteristics curve (AUC). We simulated treatment initiation over
follow up $A$ based on the logistic model $\operatorname{Pr}[A=1\mid
X]=\operatorname{expit}(0.5-2X_{1}+3X_{1}^{2}+2X_{2}-X_{3})$, where $X$ is now
a vector of predictors drawn from a 3-dimensional multivariate normal with
mean vector $\mu=(0.2,0,0.5)$ and covariance matrix
$\Sigma=\text{diag}(0.2,0.2,0.2)$. This resulted in expected treatment
initiation over follow up of 55%. We also simulated a binary outcome from a
Bernoulli distribution with mean
$\operatorname{expit}(0.2+3X_{1}-2X_{1}^{2}+2X_{2}+X_{3}-2A)$, implying an
average probability of the outcome of 66% among the untreated and 32% among
treated. Again, we set the total sample size to 1000 and randomly split the
data in a 1:1 ratio into a training and a test set.

$\displaystyle X$ $\displaystyle\sim\text{MVN}(\mu,\Sigma)$ $\displaystyle A$
$\displaystyle\sim\text{Bernoulli}\left\\{\text{expit}\left(0.5-2X_{1}+3X_{1}^{2}+2X_{2}-X_{3}\right)\right\\}$
$\displaystyle Y$
$\displaystyle\sim\text{Bernoulli}\left\\{\text{expit}\left(0.2+3X_{1}-2X_{1}^{2}+2X_{2}+X_{3}-2A\right)\right\\}$

Our prediction model was a main effects logistic regression model fit in the
training data, i.e.
$\mu_{\beta}\left(X^{*}\right)=\operatorname{expit}(\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\beta_{3}X_{3})$.
This model was misspecified with respect to the true data generating process.
We assessed the counterfactual performance of the model in an untreated
population using the AUC and the MSE, which for a binary outcome is equivalent
to the Brier score [30]. In general, positing a parametric model for
$h_{0}(X)=\mathrm{E}[(Y-g\left(X^{*}\right))^{2}\mid X,A=0]$ may be difficult
as the outcome is a squared difference. For binary outcomes, however,
expanding the square shows that to estimate $h_{0}(X)$ it is enough to
estimate $\operatorname{Pr}[Y=1\mid X,A=0]$, which is the approach we used. To
determine the effect of the specification of nuisance models $e_{a}(X)$ and
$h_{a}(X)$ on performance estimates, we compared four estimators of AUC and
MSE (${\psi}_{Na\ddot{i}ve}(0)$, ${\psi}_{IPW}(0)$, ${\psi}_{CL}(0)$, and
${\psi}_{DR}(0)$) using different combinations of correctly specified and
misspecified models for $e_{a}(X)$ and $h_{a}(X)$:

  1. 1.

Correct $e_{a}(X)$ \- main effects logistic regression model with linear and
quadratic terms.

  2. 2.

Misspecified $e_{a}(X)$ \- main effects logistic regression model with linear
terms only terms.

  3. 3.

Correct $h_{a}(X)$ \- main effects logistic regression model with linear and
quadratic terms.

  4. 4.

Misspecified $h_{a}(X)$ \- main effects logistic regression model with linear
terms only terms.

Finally, we also considered using more flexible estimation techniques for
nuisance terms $e_{a}(X)$ and $h_{a}(X)$. Specifically, we fit generalized
additive models for both using the mgcv package in $\mathrm{R}$ entering all
covariates as splines using the default options in the gam function.

Table 2: Simulation results for experiment 2 comparing the performance of the
naïve, conditional loss (CL), inverse-probability weighting (IPW), and doubly
robust (DR) estimators of the MSE and AUC under correct and misspecified
nuisance models.

| MSE | AUC
---|---|---
Estimator | Mean | $\sqrt{n}\times\text{SD}$ | $\sqrt{n}\times\text{Bias}$ | Percent | Mean | $\sqrt{n}\times\text{SD}$ | $\sqrt{n}\times\text{Bias}$ | Percent
Naïve | 0.207 | 0.176 | -0.140 | -2.1 | 0.742 | 0.491 | -1.335 | -5.4
Correct
CL | 0.212 | 0.333 | 0.015 | 0.2 | 0.783 | 0.767 | -0.045 |  1 -0.2
IPW | 0.212 | 0.517 | 0.011 | 0.2 | 0.782 | 1.258 | -0.062 |  1 -0.3
DR | 0.211 | 0.454 | 0.000 | 0.0 | 0.783 | 1.192 | -0.028 | -0.1
$e_{a}(X)$ misspecified
CL | 0.212 | 0.333 | 0.015 | 0.2 | 0.783 | 0.767 | -0.045 | -0.2
IPW | 0.221 | 0.358 | 0.316 | 4.7 | 0.762 | 0.876 | -0.699 | -2.8
DR | 0.212 | 0.349 | 0.016 | 0.2 | 0.782 | 0.841 | -0.066 | -0.3
$h_{a}(X)$ misspecified
CL | 0.217 | 0.356 | 0.194 | 2.9 | 0.777 | 0.803 | -0.224 | -0.9
IPW | 0.212 | 0.517 | 0.011 | 0.2 | 0.782 | 1.258 | -0.062 | -0.3
DR | 0.211 | 0.625 | 0.001 | 0.0 | 0.783 | 1.317 | -0.024 | -0.1
Both misspecified
CL gam | 0.213 | 0.348 | 0.052 | 0.8 | 0.782 | 0.800 | -0.063 | -0.3
IPW gam | 0.214 | 0.422 | 0.075 | 1.1 | 0.778 | 1.032 | -0.181 | -0.7
DR gam | 0.211 | 0.403 | 0.010 | 0.1 | 0.784 | 0.966 | -0.021 | -0.1
Truth | 0.211 |  |  |  | 0.784 |  |  |

  * •

Average of estimates, estimated bias, estimated standard deviation (SD), and
estimated relative bias for the naïve empirical, weighting (IPW), conditional
loss (CL), and doubly robust (DR) estimators. $\sqrt{n}$ is the number of
observations in the test set. Here, $h_{a}(X)$ is a model for
$\operatorname{Pr}[Y=1\mid X,A=a]$ and $e_{a}(X)$ denotes a model for
$\Pr[A=a|X]$. Relative bias is calculated as
$(\text{estimator}-\text{truth})/\text{truth}$. Correct and Misspecified refer
to the nuisance models, $e_{a}(X)$ or $h_{a}(X)$ or both. In the final rows,
gam indicates that a generalized additive model was used to estimate nuisance
models. Results were averaged over 10,000 simulations.

Table 2 shows the results from experiment 2. As in the previous experiment,
the naïve empirical estimators of the AUC and MSE were biased relative to the
true counterfactual values with a relative bias of 2.1% and $-5.4\%$
respectively. When all models were correctly specified, the weighting,
conditional loss, and doubly robust estimators were all unbiased (absolute
relative bias between 0.2% to 0.3%). When $e_{a}(X)$ was misspecified, the
weighting estimator was biased (relative bias of 4.7% and -2.8%) but the
conditional loss and doubly robust estimator were unbiased (absolute relative
bias of 0.2% to 0.3%). Under misspecification of $h_{a}(X)$, the conditional
loss estimator was biased (relative bias of 2.9% and -0.9%), but the weighting
estimator and the doubly robust estimator were unbiased (absolute relative
bias of 0.0% to 0.3%). When both models $e_{a}(X)$ and $h_{a}(X)$ were
misspecified all estimators, including the doubly robust estimator, were
biased. Finally, when a more flexible generalized additive model was used to
estimate both $e_{a}(X)$ and $h_{a}(X)$, the doubly robust estimator was
unbiased (absolute relative bias of 0.1%). Across all scenarios, the weighting
estimator generally had the largest standard errors and the conditional loss
estimator had the smallest standard errors.

##  9 Application to prediction of statin-naïve risk

We apply the proposed methods to evaluate the performance of two
counterfactual prediction models targeting the statin-naïve risk of
cardiovascular disease: that is the risk in the same population if, contrary
to fact, statins had been withheld. We compare one model that was explicitly
tailored for the counterfactual estimand of interest and a second that was
not.

###  9.1 Study design and data

The Multi-Ethnic Study on Atherosclerosis (MESA) study is a population-based
sample of 6,814 men and women aged 45 to 84 drawn from six communities
(Baltimore; Chicago; Forsyth County, North Carolina; Los Angeles; New York;
and St. Paul, Minnesota) in the United States between 2000 and 2002. The
design, sampling procedures, and collection methods of the study have been
described previously [31]. Study teams conducted five examination visits
between 2000 and 2011 in 18 to 24 month intervals focused on the prevalence,
correlates, and progression of subclinical cardiovascular disease. These
examinations included assessments of lipid-lowering medication use (primarily
statins), as well as assessments of cardiovascular risk factors such as
systolic blood pressure, serum cholesterol, cigarette smoking, height, weight,
and diabetes.

In a previous analysis, we used data from the MESA study to emulate a trial
comparing continuous statin use versus no statins and benchmarked our results
against those from published randomized trials. To construct a model of the
statin-naïve risk, we then emulated a single arm trial in which no one started
statins over a 10-year follow up period. To determine trial eligibility, we
followed the AHA guidelines [32] on statin use which stipulate that patients
aged 40 to 75 with serum LDL cholesterol levels between 70 mg/dL and 190 mg/dL
and no history of cardiovascular disease should initiate statins if their
(statin-naïve) risk exceeds 7.5%. Therefore, we considered MESA participants
who completed the baseline examination, had no previous history of statin use,
no history of cardiovascular disease, and who met the criteria described in
the guidelines (excluding the risk threshold) as eligible to participate in
the trial. The primary endpoint was time to atherosclerotic cardiovascular
disease (ASCVD), defined as nonfatal myocardial infarction, coronary heart
disease death, or ischemic stroke.

Follow up began at the second examination cycle to enable a “wash out” period
for statin use and to ensure adequate pre-treatment covariates to control
confounding (some were taken from the first cycle and others from the second).
In the original analysis, we constructed a sequence of nested trials starting
at each exam, however here for simplicity we limited our attention to the
first trial. We used the questionnaire in examinations three through five to
determine statin initiation over the follow up period. Because the exact
timing of statin initiation was not known with precision, we estimated it by
drawing a random month between the current and previous examinations
(sensitivity analyses conducted in original study).

Of the 6,814 MESA participants who completed the baseline examination, 4,149
met the eligibility criteria for our trial emulation. There were 288 ASCVD
events and 190 non-ASCVD deaths. For simplicity, here we dropped those lost to
follow up and who first have competing events although in practice both can be
accommodated in our framework for evaluating the performance of a
counterfactual prediction model. For model training and evaluation, we further
randomly split the dataset into training and test sets of equal size.

###  9.2 Model estimation and performance

We compared two prediction models: one that was explicitly tailored to the
statin-naïve risk and a second that was not. Both models used the same
regression specification with main effects of baseline predictors commonly
used in cardiovascular risk prediction: age, sex, smoking status, diabetes
history, systolic blood pressure, anti-hypertensive medication use and total
and HDL serum cholesterol levels.

We tailored the first model for the statin-naïve risk using inverse
probability of censoring weights. In the emulated single arm trial, statin
initiation can be viewed as “non-adherence” which can be adjusted for by
inverse probability weighting, therefore we censored participants when they
initiated statins. To calculate the stablized weights, we estimated two
logistic regression models: one for the probability of remaining untreated
given past covariate history (denominator model) and one for probability of
remaining untreated given the selected baseline predictors (numerator model).
The list of covariates in the weight models are given in section F.1 of the
Appendix. To create a prediction model for the statin-naïve risk, we used the
estimated weights to fit a weighted logistic regression model conditional on
the baseline predictors of interest.

For comparison, we fit a second traditional (factual) prediction model by
regressing the observed ASCVD event indicator on the same set of baseline
predictors, but ignoring treatment initiation over the follow up period. This
approach targets the “natural course” risk (i.e., the risk under the statin
initiation policies that prevailed at the time of the study) rather than the
statin-naïve risk. We fit the model using standard logistic regression based
on maximum likelihood.

To assess the performance of the models, we estimated the naïve and
counterfactual MSE in the test set. For the latter we used the conditional
loss, inverse probability weighting, and doubly robust estimators of the MSE.
Models for the initiation of treatment $e_{a}(X)$ and for the conditional loss
$h_{a}(X)$ were implemented as main effects logistic regression models. As in
the simulation example, to estimate the conditional loss it is sufficient to
model the probability of the outcome alone. To quantify uncertainty, we used
the non-parametric bootstrap with 1000 bootstrap replicates.

Table 3: Estimated MSE in a statin-naïve population for two prediction models
using emulated trial data from MESA.

| MSE | AUC
---|---|---
Estimator | Logit | Weighted Logit | Logit | Weighted Logit
Naïve | 0.069 | 0.072 | 0.710 | 0.708
| (0.003) | (0.003) | (0.013) | (0.014)
CL | 0.086 | 0.085 | 0.719 | 0.727
| (0.005) | (0.004) | (0.015) | (0.015)
IPW | 0.109 | 0.099 | 0.753 | 0.778
| (0.013) | (0.009) | (0.025) | (0.029)
DR | 0.090 | 0.087 | 0.740 | 0.751
| (0.006) | (0.005) | (0.023) | (0.023)

  * •

The columns refer to the posited prediction model: Logit is an (unweighted)
logistic regression model and weighted logit is a logistic regression model
with inverse probability weights for remaining statin-free. The rows are the
model performance estimates of the MSE and AUC. Naïve is the empirical
estimator using factual outcomes ($\widehat{\psi}_{Na\ddot{i}ve}(0)$), CL is
the conditional loss estimator ($\widehat{\psi}_{CL}(0)$), IPW is the inverse
probability weighting estimator ($\widehat{\psi}_{IPW}(0)$), DR is the doubly-
robust estimator ($\widehat{\psi}_{DR}(0)$). Standard error estimates are
shown in parentheses obtained via 1000 bootstrap replicates.

###  9.3 Results

Table 3 shows estimates of the AUC and MSE and the associated standard errors
in a hypothetical statin-naïve population for both prediction models using the
naïve empirical, conditional loss, weighting, and doubly robust estimators.
Across both measures, the conditional loss, weighting, and doubly robust
estimators yielded estimates that were greater (30-50% for MSE, 3-10% for AUC)
than those of the naïve empirical estimator, suggesting performance of both
models in statin-naïve population is worse than in the source population. Of
the three estimators of the statin-naïve performance, the weighting estimator
had greater standard errors than the doubly robust estimator (by 10-100%) as
well as the conditional loss estimator (by 60-160%). Consistent with the first
simulation experiment, the inverse probability weighted logistic model, which
was tailored to target the statin-naïve risk, performed worse in the source
population, but had lower MSE and higher AUC values in the counterfactual
statin-naïve population. There were sizeable differences (2-3 standard errors)
in estimates of the MSE and AUC in the counterfactual statin-naïve population
across the proposed CL, IPW, and DR estimators. All three should give
equivalent results in expectation if models for the nuisance functions
$h_{a}(X^{*})$ and $e_{a}(X^{*})$ are correctly specified. Their divergence
suggests at least one (but possibly both) of the nuisance functions may be
misspecified. Finally, drawing on the results in section D in the appendix, we
estimate the counterfactual calibration of both prediction models in a statin-
naïve population using a weighted loess estimator where the weights are based
on inverse probability of remaining statin free. Figure 2 shows the results.
As expected, the weighted prediction model, which was tailored to target the
statin-naïve risk, was better calibrated than the unweighted model, which
generally underestimated the counterfactual risk.

Figure 2: Risk calibration curves for counterfactual prediction models fit
using logistic regression with and with out inverse probability weights for
statin initiation. The rug plot shows distribution of risk predictions among
those who develop ASCVD (top) and those who don’t (bottom). The black curve is
a local-regression smoothed and inverse probability weighted estimate of the
calibration curve targeting the statin-naïve risk.

##  10 Discussion

Many practical prediction problems call for the development and evaluation of
counterfactual prediction models, for example, when the treatment distribution
changes between the time training data are obtained and when the model will be
deployed, or when predictions are meant to inform decisions about treatment
initiation. In this paper, we considered the setting when the counterfactual
prediction model has to be developed using observational training data. We
described how to estimate counterfactual prediction models and the conditions
necessary to identify these models. Separately, we also discussed how to
adjust common measures of model performance to estimate the performance of
counterfactual prediction model. Importantly, results for performance measures
were valid even when the prediction model is misspecified. A key insight was
that performance measures that have not been tailored for counterfactual
prediction will be biased. We also showed that performance can be assessed
even if the prediction model itself is misspecified. We proposed estimators
for these performance measures based on modeling the conditional loss, the
probability of treatment, and a doubly robust estimator that can be used with
data-adaptive estimators of either nuisance function.

We focused on measures of performance under a particular treatment policy.
However, prediction models may instead target the estimation of conditional
treatment effects, that is, comparisons between treatment regimes such as
$\tau(X^{*})=\operatorname{E}[Y^{1}-Y^{0}\mid X^{*}]$. In some cases, effects
may be easier to communicate to end users or may be desirable to evaluate
benefits versus harms of treatment initiation [12]. However, absolute means
and risks are common outputs of existing prediction models. Several authors
have proposed model performance measures for conditional average treatment
effects which are identified under similar assumptions to our own [33, 34, 35,
36, 37]. From an estimation standpoint, methods for targeting the conditional
average treatment effect and their performance have to balance the estimation
of the conditional risk function for the outcome under different treatment
levels with the estimation of the treatment effect function, with optimality
depending on the relative smoothness of these functions [38].

Throughout, we did not assume that the covariates needed to satisfy the
exchangeability assumption were the same covariates used in the prediction
model. This is an important aspect of our work because in practice predictors
are often chosen based on their availability in a clinical setting rather than
what would be optimal from a causal (or even predictive) perspective [19, 39].
Moreover, our methods reflect the fact that confounding is a problem that
needs to be adjusted for in the setting where the model is developed, whereas
predictions often need to be optimized for the setting where the model will be
applied (although as discussed below allowing the covariate distribution to
differ between settings can introduce further subtleties).

One limitation of our approach is that it require that a set of covariates
sufficient to ensure exchanageability can be identified at the time of
training the model. Violations of this exchangeability condition can be
examined in sensitivity analyses [40, 41], for example, to explore how
violations of this assumption might affect estimates of model performance.
Further work may also examine identification of counterfactual prediction
models and their performance under alternative identifiability conditions.

In this work, we assumed that the distribution of predictors is the same in
the training and deployment setting. In many cases, however, the covariate
distributions are also likely to differ between settings [42, 43]. Like
differences in treatment initiation, differences in covariate distributions
may cause the performance of the prediction model to degrade, particularly
when the model is misspecified. Methods for transporting prediction models
from source to target populations which relate to our own have previously been
proposed [43, 44, 45, 28, 46, 29] as have methods for transporting conditional
average effects from trials to target populations with different covariate
distributions [47, 48, 49, 50]. In future work, our results could be extended
to allow for both differences in the distribution of treatment $f(A|X)$ and
the covariates $f(X)$ between training and deployment settings.

## References

  * [1] Nan Geloven et al.  “Prediction Meets Causal Inference: The Role of Treatment in Clinical Prediction Models”  In _Eur J Epidemiol_ 35.7, 2020, pp. 619–630  DOI: 10.1007/s10654-020-00636-1
  * [2] Samuel G. Finlayson et al.  “The Clinician and Dataset Shift in Artificial Intelligence”  In _N Engl J Med_ 385.3 Massachusetts Medical Society, 2021, pp. 283–286  DOI: 10.1056/NEJMc2104626
  * [3] Adarsh Subbaswamy and Suchi Saria  “From Development to Deployment: Dataset Shift, Causality, and Shift-Stable Models in Health AI”  In _Biostatistics_ 21.2, 2020, pp. 345–352  DOI: 10.1093/biostatistics/kxz041
  * [4] Romin Pajouheshnia et al.  “Accounting for Treatment Use When Validating a Prognostic Model: A Simulation Study”  In _BMC Med Res Methodol_ 17.1, 2017, pp. 103  DOI: 10.1186/s12874-017-0375-8
  * [5] Rolf H. H. Groenwold et al.  “Explicit Inclusion of Treatment in Prognostic Modeling Was Recommended in Observational and Randomized Settings”  In _J Clin Epidemiol_ 78, 2016, pp. 90–100  DOI: 10.1016/j.jclinepi.2016.03.017
  * [6] Amanda Coston, Alan Mishler, Edward H. Kennedy and Alexandra Chouldechova  “Counterfactual Risk Assessments, Evaluation, and Fairness”, 2020  DOI: 10.48550/arXiv.1909.00066
  * [7] Barbra A. Dickerman et al.  “Predicting Counterfactual Risks under Hypothetical Treatment Strategies: An Application to HIV”  In _Eur J Epidemiol_ 37.4, 2022, pp. 367–376  DOI: 10.1007/s10654-022-00855-8
  * [8] Lijing Lin et al.  “A Scoping Review of Causal Methods Enabling Predictions under Hypothetical Interventions”  In _Diagn Progn Res_ 5.1, 2021, pp. 3  DOI: 10.1186/s41512-021-00092-9
  * [9] Peter Schulam and Suchi Saria  “Reliable Decision Support Using Counterfactual Models”  In _Advances in Neural Information Processing Systems_ 30 Curran Associates, Inc., 2017  URL: https://proceedings.neurips.cc/paper/2017/hash/299a23a2291e2126b91d54f3601ec162-Abstract.html
  * [10] Jeroen Hoogland et al.  “A Tutorial on Individualized Treatment Effect Prediction from Randomized Trials with a Binary Endpoint”  In _Statistics in Medicine_ 40.26, 2021, pp. 5961–5981  DOI: 10.1002/sim.9154
  * [11] David M. Kent, Ewout Steyerberg and David Klaveren  “Personalized Evidence Based Medicine: Predictive Approaches to Heterogeneous Treatment Effects”  In _BMJ_ 363 British Medical Journal Publishing Group, 2018, pp. k4245  DOI: 10.1136/bmj.k4245
  * [12] David M. Kent et al.  “The Predictive Approaches to Treatment Effect Heterogeneity (PATH) Statement”  In _Ann Intern Med_ 172.1 American College of Physicians, 2020, pp. 35–45  DOI: 10.7326/M18-3667
  * [13] Paul P. Glasziou and Les M. Irwig  “An Evidence Based Approach to Individualising Treatment”  In _BMJ_ 311.7016 British Medical Journal Publishing Group, 1995, pp. 1356–1359  DOI: 10.1136/bmj.311.7016.1356
  * [14] Issa J Dahabreh, Rodney Hayward and David M Kent  “Using Group Data to Treat Individuals: Understanding Heterogeneous Treatment Effects in the Age of Precision Medicine and Patient-Centred Evidence”  In _International Journal of Epidemiology_ 45.6, 2016, pp. 2184–2193  DOI: 10.1093/ije/dyw125
  * [15] Barbra A. Dickerman and Miguel A. Hernán  “Counterfactual Prediction Is Not Only for Causal Inference”  In _Eur J Epidemiol_ 35.7, 2020, pp. 615–617  DOI: 10.1007/s10654-020-00659-8
  * [16] Miguel A Hernán and James M Robins  “Causal Inference: What If”  Boca Raton: Chapman & Hall/CRC, 2020
  * [17] Frank E. Harrell, Kerry L. Lee and Daniel B. Mark  “Multivariable Prognostic Models: Issues in Developing Models, Evaluating Assumptions and Adequacy, and Measuring and Reducing Errors”  In _Statistics in Medicine_ 15.4, 1996, pp. 361–387  DOI: 10.1002/(SICI)1097-0258(19960229)15:4¡361::AID-SIM168¿3.0.CO;2-4
  * [18] Douglas G. Altman and Patrick Royston  “What Do We Mean by Validating a Prognostic Model?”  In _Statistics in Medicine_ 19.4, 2000, pp. 453–473  DOI: 10.1002/(SICI)1097-0258(20000229)19:4¡453::AID-SIM350¿3.0.CO;2-5
  * [19] Ewout W. Steyerberg  “Clinical Prediction Models: A Practical Approach to Development, Validation, and Updating”, Statistics for Biology and Health  Cham: Springer International Publishing, 2019  DOI: 10.1007/978-3-030-16399-0
  * [20] James Robins  “A New Approach to Causal Inference in Mortality Studies with a Sustained Exposure Period—Application to Control of the Healthy Worker Survivor Effect”  In _Mathematical Modelling_ 7.9, 1986, pp. 1393–1512  DOI: 10.1016/0270-0255(86)90088-6
  * [21] James Robins  “A Graphical Approach to the Identification and Estimation of Causal Parameters in Mortality Studies with Sustained Exposure Periods”  In _Journal of Chronic Diseases_ 40, 1987, pp. 139S–161S  DOI: 10.1016/S0021-9681(87)80018-8
  * [22] Donald B. Rubin  “Randomization Analysis of Experimental Data: The Fisher Randomization Test Comment”  In _Journal of the American Statistical Association_ 75.371 [American Statistical Association, Taylor & Francis, Ltd.], 1980, pp. 591–593  DOI: 10.2307/2287653
  * [23] Matthew Sperrin et al.  “Using Marginal Structural Models to Adjust for Treatment Drop-in When Developing Clinical Prediction Models”  In _Statistics in Medicine_ 37.28, 2018, pp. 4142–4154  DOI: 10.1002/sim.7913
  * [24] James M. Robins, Steven D. Mark and Whitney K. Newey  “Estimating Exposure Effects by Modelling the Expectation of Exposure Conditional on Confounders”  In _Biometrics_ 48.2 [Wiley, International Biometric Society], 1992, pp. 479–495  DOI: 10.2307/2532304
  * [25] Leonard E. Braitman and Paul R. Rosenbaum  “Rare Outcomes, Common Treatments: Analytic Strategies Using Propensity Scores”  In _Ann Intern Med_ 137.8 American College of Physicians, 2002, pp. 693–695  DOI: 10.7326/0003-4819-137-8-200210150-00015
  * [26] Victor Chernozhukov et al.  “Double/Debiased Machine Learning for Treatment and Structural Parameters”  In _Econom J_ 21.1 Oxford Academic, 2018, pp. C1–C68  DOI: 10.1111/ectj.12097
  * [27] James Robins, Lingling Li, Eric Tchetgen and Aad Vaart  “Higher Order Influence Functions and Minimax Estimation of Nonlinear Functionals”  In _Probability and Statistics: Essays in Honor of David A. Freedman_ 2 Institute of Mathematical Statistics, 2008, pp. 335–422  DOI: 10.1214/193940307000000527
  * [28] Jon A Steingrimsson, Constantine Gatsonis, Bing Li and Issa J Dahabreh  “Transporting a Prediction Model for Use in a New Target Population”  In _American Journal of Epidemiology_ 192.2, 2023, pp. 296–304  DOI: 10.1093/aje/kwac128
  * [29] Samantha Morrison et al.  “Robust Estimation of Loss-Based Measures of Model Performance under Covariate Shift”, 2022  arXiv: http://arxiv.org/abs/2210.01980
  * [30] Glenn W. Brier  “VERIFICATION OF FORECASTS EXPRESSED IN TERMS OF PROBABILITY”  In _Mon. Wea. Rev._ 78.1 American Meteorological Society, 1950, pp. 1–3  DOI: 10.1175/1520-0493(1950)078¡0001:VOFEIT¿2.0.CO;2
  * [31] Diane E. Bild et al.  “Multi-Ethnic Study of Atherosclerosis: Objectives and Design”  In _American Journal of Epidemiology_ 156.9, 2002, pp. 871–881  DOI: 10.1093/aje/kwf113
  * [32] Grundy Scott M. et al.  “2018 AHA/ACC/AACVPR/AAPA/ABC/ACPM/ADA/AGS /APhA/ASPC/NLA/PCNA Guideline on the Management of Blood Cholesterol: A Report of the American College of Cardiology/American Heart Association Task Force on Clinical Practice Guidelines”  In _Circulation_ 139.25 American Heart Association, 2019, pp. e1082–e1143  DOI: 10.1161/CIR.0000000000000625
  * [33] Alejandro Schuler, Michael Baiocchi, Robert Tibshirani and Nigam Shah  “A Comparison of Methods for Model Selection When Estimating Individual Treatment Effects”, 2018  arXiv: http://arxiv.org/abs/1804.05146
  * [34] Craig A Rolling and Yuhong Yang  “Model Selection for Estimating Treatment Effects”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 76.4 Wiley Online Library, 2014, pp. 749–769
  * [35] Yizhe Xu and Steve Yadlowsky  “Calibration Error for Heterogeneous Treatment Effects”, 2022  arXiv: http://arxiv.org/abs/2203.13364
  * [36] Mark J Laan and James M Robins  “Unified Methods for Censored Longitudinal Data and Causality”  Springer, 2003
  * [37] Ahmed Alaa and Mihaela Van Der Schaar  “Validating Causal Inference Models via Influence Functions”  In _Proceedings of the 36th International Conference on Machine Learning_ PMLR, 2019, pp. 191–201  URL: https://proceedings.mlr.press/v97/alaa19a.html
  * [38] Edward H. Kennedy  “Towards Optimal Doubly Robust Estimation of Heterogeneous Causal Effects”, 2022  arXiv: http://arxiv.org/abs/2004.14497
  * [39] Amanda Coston, Edward H. Kennedy and Alexandra Chouldechova  “Counterfactual Predictions under Runtime Confounding”, 2021  arXiv: http://arxiv.org/abs/2006.16916
  * [40] James M. Robins, Andrea Rotnitzky and Daniel O. Scharfstein  “Sensitivity Analysis for Selection Bias and Unmeasured Confounding in Missing Data and Causal Inference Models”  In _Statistical Models in Epidemiology, the Environment, and Clinical Trials_ , The IMA Volumes in Mathematics and Its Applications  New York, NY: Springer, 2000, pp. 1–94  DOI: 10.1007/978-1-4612-1284-3˙1
  * [41] Jon A. Steingrimsson, Sarah E. Robertson and Issa J. Dahabreh  “Sensitivity Analysis for Studies Transporting Prediction Models”, 2023  arXiv: http://arxiv.org/abs/2306.08084
  * [42] Steffen Bickel, Michael Brückner and Tobias Scheffer  “Discriminative Learning Under Covariate Shift”  In _J. Mach. Learn. Res._ 10, 2009, pp. 2137–2155
  * [43] Masashi Sugiyama, Matthias Krauledat and Klaus-Robert Müller  “Covariate Shift Adaptation by Importance Weighted Cross Validation”  In _Journal of Machine Learning Research_ 8.35, 2007, pp. 985–1005  URL: http://jmlr.org/papers/v8/sugiyama07a.html
  * [44] Steffen Bickel, Michael Brückner and Tobias Scheffer  “Discriminative Learning for Differing Training and Test Distributions”  In _Proceedings of the 24th International Conference on Machine Learning_ , ICML ’07  New York, NY, USA: Association for Computing Machinery, 2007, pp. 81–88  DOI: 10.1145/1273496.1273507
  * [45] Masashi Sugiyama and Motoaki Kawanabe  “Machine Learning in Non-Stationary Environments: Introduction to Covariate Shift Adaptation”  MIT press, 2012
  * [46] Bing Li, Constantine Gatsonis, Issa J. Dahabreh and Jon A. Steingrimsson  “Estimating the Area under the ROC Curve When Transporting a Prediction Model to a Target Population”  In _Biometrics_ , 2022, pp. biom.13796  DOI: 10.1111/biom.13796
  * [47] Megha L Mehrotra et al.  “Transporting Subgroup Analyses of Randomized Controlled Trials for Planning Implementation of New Interventions”  In _American Journal of Epidemiology_ 190.8, 2021, pp. 1671–1680  DOI: 10.1093/aje/kwab045
  * [48] Marissa J. Seamans et al.  “Generalizability of Subgroup Effects”  In _Epidemiology_ 32.3, 2021, pp. 389–392  DOI: 10.1097/EDE.0000000000001329
  * [49] Sarah E. Robertson et al.  “Estimating Subgroup Effects in Generalizability and Transportability Analyses”, 2021  arXiv: http://arxiv.org/abs/2109.14075
  * [50] Sarah E. Robertson, Jon A. Steingrimsson and Issa J. Dahabreh  “Regression-Based Estimation of Heterogeneous Treatment Effects When Extending Inferences from a Randomized Trial to a Target Population”  In _European Journal of Epidemiology_ 38.2, 2023, pp. 123–133  DOI: 10.1007/s10654-022-00901-5

##  Appendix A Time-fixed treatment initiation

###  A.1 Tailoring models for counterfactual predictions

Our goal is to build a model that targets the expected potential outcome under
a hypothetical intervention, e.g. the parametric model

$\operatorname{E}[Y^{a}\mid X^{*}]=\mu_{\beta}(X^{*}).$

However, we do not observe $Y^{a}$ for all individuals. Nonetheless, as the
following theorem shows, targets like $\operatorname{E}[Y^{a}\mid X^{*}]$ are
identifiable from the observed data $(X,A,Y)$ under certain assumptions.

###### Theorem 1

Under conditions A1-A3 in section 4, $\operatorname{E}[Y^{a}\mid X^{*}]$ is
identified by the observed data functionals

$\operatorname{E}[Y^{a}\mid X^{*}]=\operatorname{E}[\operatorname{E}[Y\mid
X,A=a,D_{train}=1]\mid X^{*},D_{train}=1]$ (A1)

and

$\operatorname{E}[Y^{a}\mid
X^{*}]=\operatorname{E}\left[\frac{I(A=a)}{\Pr[A=a\mid X,D_{train}=1]}Y\mid
X^{*},D_{train}=1\right]$ (A2)

in which case we can build a model for $\operatorname{E}[Y^{a}\mid X^{*}]$ by
targeting either estimand in the training dataset.

Proof. For the first representation we have

$\displaystyle\operatorname{E}[Y^{a}\mid X^{*}]$
$\displaystyle=\operatorname{E}[Y^{a}\mid X^{*},D_{train}=1]$
$\displaystyle=\operatorname{E}(\operatorname{E}[Y^{a}\mid X,D_{train}=1]\mid
X^{*},D_{train}=1)$ $\displaystyle=\operatorname{E}(\operatorname{E}[Y^{a}\mid
X,A=a,D_{train}=1]\mid X^{*},D_{train}=1)$
$\displaystyle=\operatorname{E}(\operatorname{E}[Y\mid X,A=a,D_{train}=1]\mid
X^{*},D_{train}=1)$

where the first line follows from the random sampling of the training set, the
second from the law of iterated expectations, the third from the
exchangeability condition A1, and the fourth from the consistency condition
A2. Recall that $X^{*}$ is a subset of $X$. For the second representation, we
show that it is equivalent to the first

$\displaystyle\operatorname{E}[Y^{a}\mid X^{*}]$
$\displaystyle=\operatorname{E}(\operatorname{E}[Y\mid X,A=a,D_{train}=1]\mid
X^{*},D_{train}=1)$
$\displaystyle=\operatorname{E}\left(\operatorname{E}\left[\frac{I(A=a)}{\Pr[A=a\mid
X,D_{train}=1]}Y\mid X,D_{train}=1\right]\mid X^{*},D_{train}=1\right)$
$\displaystyle=\operatorname{E}\left[\frac{I(A=a)}{\Pr[A=a\mid
X,D_{train}=1]}Y\mid X^{*},D_{train}=1\right]$

where the second line follows from the definition of conditional expectation,
and the last reverses the law of iterated expectations.

###  A.2 Identification of general loss functions

Many common model performance measures, such as the mean squared error, Brier
score, and absolute error are special cases of a generic loss function
$L(Y,\mu_{\widehat{\beta}}(X^{*}))$. To assess the performance of
counterfactual predictions, we would like to estimate

$\psi_{\widehat{\beta}}(a)=\operatorname{E}[L(Y^{a},\mu_{\widehat{\beta}}(X^{*}))]$

where $Y^{a}$ is not observed for all individuals. The following theorem
states that, under the conditions of section 4, $\psi_{\widehat{\beta}}(a)$ is
identified using the observed data alone.

###### Theorem 2

Under conditions A1-A3 in section 4 and the time-fixed setup described in
section 2, the expected loss is identified by the functionals

$\psi_{\widehat{\beta}}(a)=\operatorname{E}\left(\operatorname{E}[L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid
X,A=a,D_{test}=1]\mid D_{test}=1\right)$ (A3)

and

$\psi_{\widehat{\beta}}(a)=\operatorname{E}\left[\frac{I(A=a)}{\Pr[A=a\mid
X,D_{test}=1]}L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid D_{test}=1\right]$
(A4)

in the test set for generic counterfactual loss function
$L(Y^{a},\mu_{\widehat{\beta}})$.

Proof. For the first representation we have

$\displaystyle\psi_{\widehat{\beta}}(a)$
$\displaystyle=\operatorname{E}[L\\{Y^{a},\mu_{\widehat{\beta}}(X^{*})\\}]$
$\displaystyle=\operatorname{E}[L\\{Y^{a},\mu_{\widehat{\beta}}(X^{*})\\}\mid
D_{test}=1]$
$\displaystyle=\operatorname{E}(\operatorname{E}[L\\{Y^{a},\mu_{\widehat{\beta}}(X^{*})\\}\mid
X,D_{test}=1]\mid D_{test}=1)$
$\displaystyle=\operatorname{E}(\operatorname{E}[L\\{Y^{a},\mu_{\widehat{\beta}}(X^{*})\\}\mid
X,A=a,D_{test}=1]\mid D_{test}=1)$
$\displaystyle=\operatorname{E}(\operatorname{E}[L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid
X,A=a,D_{test}=1]\mid D_{test}=1)$

where the first line follows from the definition of
$\psi_{\widehat{\beta}}(a)$, the second from random sampling of the test set,
the third from the law of iterated expectations, the fourth from the
exchangeability condition A1, and the fifth from the consistency condition A2.
Recall that $X^{*}$ is a subset of $X$. For the second representation, we show
that it is equivalent to the first

$\displaystyle\psi_{\widehat{\beta}}(a)$
$\displaystyle=\operatorname{E}(\operatorname{E}[L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid
X,A=a,D_{test}=1]\mid D_{test}=1)$
$\displaystyle=\operatorname{E}\left(\operatorname{E}\left[\frac{I(A=a)}{\Pr[A=a\mid
X,D_{test}=1]}L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid
X,D_{test}=1\right]\mid D_{test}=1\right)$
$\displaystyle=\operatorname{E}\left[\frac{I(A=a)}{\Pr[A=a\mid
X,D_{test}=1]}L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid D_{test}=1\right]$

where the second line follows from the definition of conditional expectation,
and the last reverses the law of iterated expectations.

###  A.3 Plug-in estimation

Using sample analogs for the identified expressions A3 and A4, we obtain two
plug-in estimators for the expected loss for a generalized loss function

$\widehat{\psi}_{CL}(a)=\frac{1}{n_{test}}\sum_{i=1}^{n}I(D_{test,i}=1)\widehat{h}_{a}(X_{i})$

and

$\widehat{\psi}_{IPW}(a)=\frac{1}{n_{test}}\sum_{i=1}^{n}\frac{I(A_{i}=a,D_{test,i}=1)}{\widehat{e}_{a}(X_{i})}L\\{Y_{i},\mu_{\widehat{\beta}}(X^{*}_{i})\\}$

where $\widehat{h}_{a}(X)$ is an estimator for
$\operatorname{E}[L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid X,A=a,D_{test}=1]$
and $\widehat{e}_{a}(X)$ is an estimator for $\Pr[A=a\mid X,D_{test}=1]$.
Using the terminology in Morrison et al., we call the first plug-in estimator
the conditional loss estimator $\widehat{\psi}_{CL}(a)$ and the second the
inverse probability weighted estimator $\widehat{\psi}_{IPW}(a)$.

###  A.4 Random and dynamic regimes

Above we consider static interventions which set treatment $A$ to a particular
value $a$. We might also consider interventions which probabilistically set
$A$ based on a known density, possibly conditional on pre-treatment
covariates, e.g. $f^{*}(A\mid X)$. For instance, instead of a counterfactual
prediction if everyone or no one had been treated, we may be interested in the
prediction if 20% or 50% were treated. We term such an intervention a random
intervention to contrast it with static interventions considered previously.
Random interventions are closer to the counterfactual interventions of
interest under dataset shift which may be approximated as probabilistic
changes in the natural course of treatment due to changes in guidelines or
prescribing patterns or the wider-availability. For general counterfactual
loss function $L\\{Y^{g},\mu_{\widehat{\beta}}\\}$, the expected loss under a
random intervention is identified by the functionals

$\psi_{\widehat{\beta}}(A)=\operatorname{E}\left\\{\operatorname{E}_{f^{*}(A|X)}\left(\operatorname{E}[L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid
X,A,D_{test}=1]\mid D_{test}=1\right),D_{test}=1\right\\}$ (A5)

and

$\psi_{\widehat{\beta}}(A)=\operatorname{E}\left[\frac{f^{*}(A|X)}{f(A\mid
X,D_{test}=1)}L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid D_{test}=1\right]$
(A6)

in the test set under the time-fixed setup described in section 2. The primary
difference between these expressions and the ones in section A.1. is that the
expectation is taken with respect to the intervention density.

##  Appendix B Time-varying treatment initiation

###  B.1 Set up

Here we extend the set up of section 2 in the case that treatment initiation
is time-varying over the follow up period. We now observe $n$ i.i.d.
longitudinal samples $\\{O_{i}\\}_{i=1}^{n}$ from a source population. For
each observation, let

$O=(\overline{X}_{K},\overline{A}_{K},Y_{K+1})$

where $X_{k}$ is a vector of covariates, $A_{k}$ is an indicator of treatment,
and $Y_{k}$ is an event indicator all measured at time $k$, where
$k\in\\{0,1,\ldots,K\\}$. Overbars denote the full history of a variable, such
that $\overline{X}_{k}=(X_{0},\dots,X_{k})$. Again we assume the goal is to
build a prediction model for end of follow up outcome $Y_{K+1}$ conditional on
baseline covariates $X^{*}$ which are now a subset of baseline covariates
$X_{0}$, i.e. $X^{*}\subset X_{0}$. An example DAG for a two time point
process is shown in Figure A1.

We would like to assess the performance of the model in a counterfactual
version of the source population in which a new treatment policy is
implemented. As previously, $Y^{a}$ is the potential outcome under an
intervention which sets treatment $A$ to $a$. For a sequence of time-varying
treatments $\overline{A}_{k}$, we further define a treatment regime as a
collection of functions
$\\{g_{k}(\overline{a}_{k-1},\overline{x}_{k}):k=0,\ldots,K\\}$ for
determining treatment assignment at each time $k$, possibly based on past
treatment and covariate history. For a hypothetical treatment regime $g$, we
would like to determine the performance of fitted model
$\mu_{\widehat{\beta}}(X^{*})$ under the new regime by estimating the expected
loss

$\psi_{\widehat{\beta}}(g)=\operatorname{E}[L\\{Y^{g},\mu_{\widehat{\beta}}(X^{*})\\}]$

for generalized loss function $L\\{Y^{g},\mu_{\widehat{\beta}}(X^{*})\\}$.

$L_{0}$$A_{0}$$L_{1}$$A_{1}$$Y_{2}$$U$$P_{0}$$P_{1}$ (a) Example two time
point directed acyclic graph for prediction.

$L_{0}$$A_{0}\mid a_{0}$$L^{a_{0}}_{1}$$A^{a_{0}}_{1}\mid
a_{1}$$Y^{a_{0},a_{1}}_{2}$$U$$P_{0}$$P^{a_{0}}_{1}$ (b) Single world
intervention graph of intervention on $A_{0}$ and $A_{1}$.

Figure A1: Example directed acyclic graph (DAG) and single world intervention
graph (SWIG) for a two time point process.

###  B.2 Identifiability conditions

We now consider modified identifiability conditions under time-varying
treatment initiation. For all $k$ from 0 to $K$, we require

  1. B1.

Exchangeability: $Y^{g}_{K+1}\perp\\!\\!\\!\perp
A_{k}\mid\overline{X}_{k},\overline{A}_{k-1}$

  2. B2.

Consistency: $Y_{K+1}=Y^{g}_{K+1}$ and
$\overline{X}_{k}=\overline{X}^{g}_{k}$ if
$\overline{A}_{k}=\overline{a}_{k}^{g}$

  3. B3.

Positivity:
$1>\Pr[A_{k}=a_{k}\mid\overline{X}_{k}=\overline{x}_{k},\overline{A}_{k-1}=\overline{a}_{k-1}]>0$

###  B.3 Identification of general loss functions

The following theorem extends Theorem 2 to the case of a time-varying
treatment.

###### Theorem 3

Under conditions B1-B3 above, the expected counterfactual loss under time-
varying regime $g$ is identified by the functionals

$\displaystyle\begin{split}\psi_{\widehat{\beta}}(g)&=\operatorname{E}_{X_{0}}\bigg{[}\operatorname{E}_{X_{1}}\bigg{\\{}\ldots\operatorname{E}_{X_{K-1}}\bigg{(}\operatorname{E}_{X_{K}}[L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid\overline{X}_{K},\overline{A}_{K}=\overline{a}^{g}_{K},D_{test}=1]\\\
&\qquad\big{|}\;\overline{X}_{K-1},\overline{A}_{K-1}=\overline{a}^{g}_{K-1},D_{test}=1\bigg{)}\ldots\;\big{|}\;X_{0},A_{0}=a^{g}_{0},D_{test}=1\bigg{\\}}\;\big{|}\;D_{test}=1\bigg{]}\end{split}$
(A7)

and

$\psi_{\widehat{\beta}}(g)=\operatorname{E}\left[\frac{I(\overline{A}_{K}=\overline{a}^{g}_{K},D_{test}=1)}{\prod_{k=0}^{K}\Pr[A_{k}=a^{g}_{k}\mid\overline{X}_{k},\overline{A}_{k-1}=\overline{a}^{g}_{k-1},D_{test}=1]}L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid
D_{test}=1\right]$ (A8)

in the test set for general loss function
$L\\{Y^{g},\mu_{\widehat{\beta}}\\}$, where the first expression is a sequence
of iterated expectations and the second is an inverse-probability weighted
expectation.

Proof. For the first representation we have

$\displaystyle\psi_{\widehat{\beta}}(g)$
$\displaystyle=\operatorname{E}[L\\{Y^{g},\mu_{\widehat{\beta}}(X^{*})\\}]$
$\displaystyle=\operatorname{E}[L\\{Y^{g},\mu_{\widehat{\beta}}(X^{*})\\}\mid
D_{test}=1]$
$\displaystyle=\operatorname{E}(\operatorname{E}[L\\{Y^{g},\mu_{\widehat{\beta}}(X^{*})\\}\mid
X_{0},D_{test}=1]\mid D_{test}=1)$
$\displaystyle=\operatorname{E}(\operatorname{E}[L\\{Y^{g},\mu_{\widehat{\beta}}(X^{*})\\}\mid
X_{0},A_{0}=a^{g}_{0},D_{test}=1]\mid D_{test}=1)$

where the first line follows from the definition of
$\psi_{\widehat{\beta}}(g)$, the second from random sampling of the test set,
the third from the law of iterated expectations, and the fourth from the
exchangeability condition B1. Arguing recursively from $k=0$ to $K$, we can
repeatedly invoke iterated expectations and exchanageability to insert
$\overline{X}_{k}$ and $\overline{A}_{k}=\overline{a}^{g}_{k}$, such that

$\displaystyle\psi_{\widehat{\beta}}(g)$
$\displaystyle=\operatorname{E}_{X_{0}}\bigg{[}\operatorname{E}_{X_{1}}\bigg{\\{}\ldots\operatorname{E}_{X_{K}}\bigg{(}E[L\\{Y^{g},\mu_{\widehat{\beta}}(X^{*})\\}\mid\overline{X}_{K},\overline{A}_{K}=\overline{a}^{g}_{K},D_{test}=1]$
$\displaystyle\qquad\big{|}\;\overline{X}_{K},\overline{A}_{K}=\overline{a}^{g}_{K},D_{test}=1\bigg{)}\ldots\;\big{|}\;X_{0},A_{0}=a^{g}_{0},D_{test}=1\bigg{\\}}\;\big{|}\;D_{test}=1\bigg{]}$
$\displaystyle=\operatorname{E}_{X_{0}}\bigg{[}\operatorname{E}_{X_{1}}\bigg{\\{}\ldots\operatorname{E}_{X_{K}}\bigg{(}\operatorname{E}[L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid\overline{X}_{K},\overline{A}_{K}=\overline{a}^{g}_{K},D_{test}=1]$
$\displaystyle\qquad\big{|}\;\overline{X}_{K},\overline{A}_{K}=\overline{a}^{g}_{K},D_{test}=1\bigg{)}\ldots\;\big{|}\;X_{0},A_{0}=a^{g}_{0},D_{test}=1\bigg{\\}}\;\big{|}\;D_{test}=1\bigg{]}$

where the last line follows by consistency condition B2. For the second
representation, note that for the inner most expectations we can proceed as
previously

$\displaystyle\operatorname{E}(\operatorname{E}[L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid\overline{X}_{K},\overline{A}_{k}=\overline{a}^{g}_{K},D_{test}=1]\mid\overline{X}_{K-1},\overline{A}_{k-1}=\overline{a}^{g}_{K-1},D_{test}=1)$
$\displaystyle=\operatorname{E}\left(\operatorname{E}\left[W_{K}L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid\overline{X}_{K},\overline{A}_{K-1},D_{test}=1\right]\mid\overline{X}_{K-1},\overline{A}_{K-1},D_{test}=1\right)$
$\displaystyle=\operatorname{E}\left(W_{K}\operatorname{E}\left[L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid\overline{X}_{K},\overline{A}_{K-1},D_{test}=1\right]\mid\overline{X}_{K-1},\overline{A}_{K-1},D_{test}=1\right)$
$\displaystyle=\operatorname{E}\left[W_{K}L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid\overline{X}_{K-1},\overline{A}_{K-1},D_{test}=1\right]$

where the second line follows from the definition of conditional expectation,
the third removes the constant fraction outside expectation, and the last
reverses the law of iterated expectations and where

$W_{K}=\frac{I(A_{K}=a_{K}^{g},D_{test}=1)}{\Pr[A_{K}=a_{K}^{g}\mid\overline{X}_{K},\overline{A}_{K-1},D_{test}=1]}$

Arguing recursively from $k=0$ to $K$, we get

$\displaystyle\psi_{\widehat{\beta}}(g)$
$\displaystyle=\operatorname{E}\left[\frac{I(\overline{A}_{K}=\overline{a}^{g}_{K},D_{test}=1)}{\prod_{k=0}^{K}\Pr[A_{k}=a^{g}_{k}\mid\overline{X}_{k},\overline{A}_{k-1}=\overline{a}^{g}_{k-1},D_{test}=1]}L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid
D_{test}=1\right]$

which is the inverse-probability weighted representation with weights equal to

$W_{k}=\frac{I(\overline{A}_{K}=\overline{a}^{g}_{K},D_{test}=1)}{\prod_{k=0}^{K}\Pr[A_{k}=a^{g}_{k}\mid\overline{X}_{k},\overline{A}_{k-1}=\overline{a}^{g}_{k-1},D_{test}=1]}.$

###  B.4 Plug-in estimation

Using sample analogs for the identified expressions A7 and A8, we obtain two
plug-in estimators for the expected counterfactual loss under a generalized
loss function

$\widehat{\psi}_{CL}(g)=\frac{1}{n_{test}}\sum_{i=1}^{n}I(D_{test,i}=1)\widehat{h}_{a_{0}}(X_{i})$

and

$\widehat{\psi}_{IPW}(g)=\frac{1}{n_{test}}\sum_{i=1}^{n}\frac{I(\overline{A}_{K,i}=\overline{a}^{g}_{K,i},D_{test,i}=1)}{\prod_{k=0}^{K}\widehat{e}_{a_{k}}(X_{i})}L\\{Y_{i},\mu_{\widehat{\beta}}(X^{*}_{i})\\}$

where $h_{t+1}=L\\{Y,\mu_{\widehat{\beta}}(X^{*}_{i})\\}$ and $h_{a_{0}}(X)$
is recursively defined for $t$ = $K,\ldots,0$

$h_{a_{t}}:(x_{t},a_{t})\operatorname{E}[h_{a_{t+1}}(X_{t+1})\mid\overline{X}_{t},\overline{A}_{t}=\overline{a}^{g}_{t}]$

$\widehat{h}_{a}(X)$ is an estimator for
$\operatorname{E}[L\\{Y,\mu_{\widehat{\beta}}(X^{*})\\}\mid X,A=a,D_{test}=1]$
and $\widehat{e}_{a_{k}}(X)$ is an estimator for
$\Pr[A_{k}=a^{g}_{k}\mid\overline{X}_{k},\overline{A}_{k-1}=\overline{a}^{g}_{k-1},D_{test}=1]$.
Note that as the number of time points (i.e. $K$) increases, the proportion in
the test set who actually follow the regime of interest, i.e. those for whom
$I(\overline{A}_{K}=\overline{a}^{g}_{K},D_{test,i}=1)=1$ may be prohibitively
small, in which case plug-in estimation may not be feasible. In this case,
additional modeling assumptions will be necessary to borrow information from
other regimes.

##  Appendix C Doubly robust estimators

###  C.1 Efficient influence function

We have shown previously that, under the identifiability conditions of section
4, the expected counterfactual loss of a generalized loss function
$L\\{Y^{a},\mu(X^{*})\\}$ is identified by the observed data functional

$\psi(a)=\operatorname{E}\left(\operatorname{E}[L\\{Y,\mu(X^{*})\\}\mid
X,A=a]\right).$

In the following theorem, we identify the efficient influence function for
$\psi(a)$ under a nonparametric model for the observed data in the test set.

###### Theorem 4

The influence function for $\psi(a)$ under a nonparametric model for the
observable data $O=(X,A,Y)$ is

$\displaystyle\chi_{P_{0}}^{1}$ $\displaystyle=\frac{I(A=a)}{\Pr[A=a\mid
X]}(L\\{Y,\mu(X^{*})\\}-\operatorname{E}[L\\{Y,\mu(X^{*})\\}\mid X,A=a])\;+$
$\displaystyle\qquad(\operatorname{E}[L\\{Y,\mu(X^{*})\\}\mid
X,A=a]-\psi(a)).$

As the influence function under a nonparametric model is always unique, it is
also the efficient influence function.

Proof. To show that $\chi_{P_{0}}^{1}$ is the efficient influence function, we
will use the well-known fact that the influence function is a solution to

$\frac{d}{dt}\psi_{P_{t}}(a)\bigg{|}_{t=0}=\operatorname{E}_{P_{0}}(\chi_{P_{0}}^{1}g_{P_{0}})$

where $g_{P_{0}}$ is the score of the obeservable data under the true law
$P_{0}$ and $P_{t}$ is a parametric submodel indexed by $t\in[0,1]$ and the
pathwise derivative of the submodel is evaluated at $t=0$ corresponding to the
true law $P_{0}$. Let $h_{a}(X)=E_{P_{0}}[L\\{Y,\mu(X^{*})\\}\mid X,A=a]$.
Beginning with the left hand side

$\displaystyle\frac{d}{dt}\psi_{P_{t}}(a)\bigg{|}_{t=0}$
$\displaystyle=\frac{d}{dt}\operatorname{E}_{P_{t}}\left(\operatorname{E}_{P_{t}}[L\\{Y,\mu(X^{*})\\}\mid
X,A=a]\right)\bigg{|}_{t=0}$ $\displaystyle=\frac{\partial}{\partial
t}\operatorname{E}_{P_{t}}\left(\operatorname{E}_{P_{0}}[L\\{Y,\mu(X^{*})\\}\mid
X,A=a]\right)\bigg{|}_{t=0}\;+$
$\displaystyle\qquad\operatorname{E}_{P_{0}}\left(\frac{\partial}{\partial
t}\operatorname{E}_{P_{t}}[L\\{Y,\mu(X^{*})\\}\mid
X,A=a]\bigg{|}_{t=0}\right)$
$\displaystyle=\operatorname{E}_{P_{0}}\left[\left\\{h_{a}(X)-\psi(a)\right\\}g_{X,A,Y}(O)\right]\;+$
$\displaystyle\qquad\operatorname{E}_{P_{0}}\left\\{\left(\frac{I(A=a)}{\Pr[A=a\mid
X]}\bigg{[}L\\{Y,\mu(X^{*})\\}-h_{a}(X)\bigg{]}\right)g_{X,A,Y}(O)\right\\}$
$\displaystyle=\operatorname{E}_{P_{0}}\left\\{\left(h_{a}(X)-\psi(a)+\frac{I(A=a)}{\Pr[A=a\mid
X]}\bigg{[}L\\{Y,\mu(X^{*})\\}-h_{a}(X)\bigg{]}\right)g_{X,A,Y}(O)\right\\}$

where the first line is the definition, the second line applies the chain
rule, the third applies definition of the score, and the last uses linearity
of expectations. Returning to original supposition, it follows that the
influence function is

$\displaystyle\chi_{P_{0}}^{1}$ $\displaystyle=\frac{I(A=a)}{\Pr[A=a\mid
X]}(L\\{Y,\mu(X^{*})\\}-\operatorname{E}[L\\{Y,\mu(X^{*})\\}\mid X,A=a])\;+$
$\displaystyle\qquad(\operatorname{E}[L\\{Y,\mu(X^{*})\\}\mid
X,A=a]-\psi(a)).$

###  C.2 One-step estimator

Given the efficient influence function above and random sampling in the test
set, the one-step estimator for $\psi(a)$ is given by

$\widehat{\psi}_{DR}(a)=\frac{1}{n_{test}}\sum_{i=1}^{n}I(D_{test,i}=1)\widehat{h}_{a}(X_{i})+\frac{I(A_{i}=a,D_{test,i}=1)}{\widehat{e}_{a}(X_{i})}\left[L\\{Y_{i},\mu(X^{*}_{i})\\}-\widehat{h}_{a}(X_{i})\right]$

###  C.3 Asymptotic properties

In previous sections, the asymptotic properties of $\widehat{\psi}_{CL}(a)$
and $\widehat{\psi}_{IPW}(a)$ follow from standard parametric theory111after
separating estimation of $\mu_{\beta}(X^{*})$ from the evaluation of
performance by random partition of test set.. However, the asymptotic
properties of $\widehat{\psi}_{DR}(a)$ are complicated by the estimation of
two nuisance functions, $\widehat{h}_{a}(X)$ and $\widehat{e}_{a}(X)$, and the
fact that, we do not immediately assume a parametric model for either. To
simplify the derivation of the large sample properties of
$\widehat{\psi}_{DR}(a)$ we begin by defining

$H\left(e_{a}^{\prime}(X),h_{a}^{\prime}(X)\right)=h_{a}^{\prime}(X)+\frac{I(A=a)}{e_{a}^{\prime}(X)}\left[L\left(Y,\mu\left(X^{*}\right)\right)-h_{a}^{\prime}(X)\right]$

for arbitrary functions $e_{a}^{\prime}(X)$, and $h_{a}^{\prime}(X)$. Here we
suppress the dependence on being in the test set for ease of exposition, but
note that the rest procedes the same if we were to limit our focus to the test
set. Note, the doubly robust estimator can be written as
$\widehat{\psi}_{DR}(a)=\frac{1}{n}\sum_{i=1}^{n}H\left(\widehat{e}_{a}\left(X_{i}\right),\widehat{h}_{a}\left(X_{i}\right)\right)$.
We define the probability limits of $\widehat{e}_{a}(X)$ and
$\widehat{h}_{a}(X)$ as $e_{a}^{*}(X)$ and $h_{a}^{*}(X)$, respectively. By
definition, when $\widehat{e}_{a}(X)$ and $\widehat{h}_{a}(X)$ are correctly
specified, the limits are $e_{a}^{*}(X)=$ $\operatorname{Pr}[A=a\mid X]$ and
$h^{*}_{a}(X)=\mathrm{E}\left[L\left(Y,\mu\left(X^{*}\right)\right)\mid
X,A=a\right]$.

To derive the asymptotic properties of $\widehat{\psi}_{DR}(a)$, we make the
following assumptions:

  1. D1.

$H(\widehat{e}_{a}(X),\widehat{h}_{a}(X))$ and its limit
$H\left(e^{*}_{a}(X),h^{*}_{a}(X)\right)$ fall in a Donsker class.

  2. D2.

$\left\|H(\widehat{e}_{a}(X),\widehat{h}_{a}(X))-H\left(e^{*}_{a}(X),h^{*}_{a}(X)\right)\right\|\stackrel{{\scriptstyle
P}}{{\longrightarrow}}0$.

  3. D3.

(Finite second moment).
$\mathrm{E}\left[H\left(e^{*}_{a}(X),h^{*}_{a}(X)\right)^{2}\right]<\infty$.

  4. D4.

(Model double robustness). At least one of the models $\widehat{e}_{a}(X)$ or
$\widehat{h}_{a}(X)$ is correctly specified. That is, at least one of
$e^{*}_{a}(X)=\operatorname{Pr}[A=a\mid X]$ or
$h^{*}_{a}(X)=\mathrm{E}\left[L\left(Y,\mu\left(X^{*}\right)\right)\mid
X,A=a\right]$ holds, but not necessarily both.

Assumption D1 is a well-known restriction on the complexity of the functionals
$\widehat{e}_{a}(X)$ and $\widehat{h}_{a}(X)$. As long as
$\widehat{e}_{a}(X),\widehat{h}_{a}(X),e^{*}_{a}(X)$, and $h^{*}_{a}(X)$ are
Donsker and all are uniformly bounded then Assumption D1 holds by the Donsker
preservation theorem. Many commonly used models such as generalized linear
models fall within the Donsker class. This requirement can be further relaxed
through sample-splitting, in which case more flexible machine learning
algorithms such as random forests, gradient boosting, or neural networks may
be used to estimate $\widehat{e}_{a}(X)$ and $\widehat{h}_{a}(X)$.

Using Assumptions D1 through D4, below we prove:

  1. 1.

(Consistency) $\widehat{\psi}_{DR}(a)\stackrel{{\scriptstyle
P}}{{\longrightarrow}}\psi(a)$.

  2. 2.

(Asymptotic distribution) $\widehat{\psi}_{DR}(a)$ has the asymptotic
representation

$\sqrt{n}\left(\widehat{\psi}_{DR}(a)-\psi(a)\right)=\sqrt{n}\left(\frac{1}{n}\sum_{i=1}^{n}H\left(e^{*}_{a}(X_{i}),h^{*}_{a}(X_{i})\right)-\mathrm{E}\left[H\left(e^{*}_{a}(X),h^{*}_{a}(X)\right)\right]\right)+Re+o_{P}(1),$

where

$Re\leq\sqrt{n}O_{P}\left(\left\|\widehat{h}_{a}(X)-\mathrm{E}\left[L\left(Y,\mu(X^{*})\right)\mid
X,A=a\right]\right\|_{2}^{2}\times\Big{\|}\widehat{e}_{a}(X)-\operatorname{Pr}[A=a\mid
X]\Big{\|}_{2}^{2}\right)$

and thus if $\widehat{h}_{a}(X)$ and $\widehat{e}_{a}(X)$ converge at combined
rate of at least $\sqrt{n}$ then

$\sqrt{n}\left(\widehat{\psi}_{DR}(a)-\psi(a)\right)\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N\left(0,\operatorname{Var}\left[H(e^{*}_{a}(X),h^{*}_{a}(X))\right]\right)$

####  C.3.1 Consistency

Using the probability limits $e_{a}^{*}(X)$ and $h_{a}^{*}(X)$ defined
previously, the double robust estimator $\widehat{\psi}_{DR}(a)$ converges in
probability to

$\widehat{\psi}_{DR}(a)\stackrel{{\scriptstyle
P}}{{\longrightarrow}}\operatorname{E}\left[h^{*}_{a}(X)+\frac{I(A=a)}{e^{*}_{a}(X)}\left(L\left(Y,\mu\left(X^{*}\right)\right)-h^{*}_{a}(X)\right)\right]$

Here we show that the right-hand side is equal to $\psi$ under assumptions D1-
D4 when either:

  1. 1.

$\widehat{e}_{a}(X)$ is correctly specified

  2. 2.

$\widehat{h}_{a}(X)$ is correctly specified

First consider the case where $\widehat{e}_{a}(X)$ is correctly specified,
that is $e^{*}_{a}(X)=\operatorname{Pr}[A=a\mid X]$, but we do not assume that
the limit $h^{*}_{a}(X)$ is equal to
$\left.\mathrm{E}\left[L\left(Y,g\left(X^{*}\right)\right)\mid
X,A=a\right]\right)$. Recall, as shown previously
$\psi=\operatorname{E}\left[\frac{I(A=a)}{\Pr[A=a\mid
X]}L(Y,\mu_{\widehat{\beta}}(X^{*}))\right]$

$\displaystyle\widehat{\psi}_{DR}(a)$ $\displaystyle\stackrel{{\scriptstyle
P}}{{\rightarrow}}\operatorname{E}\left[h^{*}_{a}(X)+\frac{I(A=a)}{e^{*}_{a}(X)}\left(L\left(Y,\mu\left(X^{*}\right)\right)-h^{*}_{a}(X)\right)\right]$
$\displaystyle=\operatorname{E}\left[h^{*}_{a}(X)-\frac{I(A=a)}{e^{*}_{a}(X)}h^{*}_{a}(X)\right]+\psi(a)$
$\displaystyle=\operatorname{E}\left[\operatorname{E}\left[h^{*}_{a}(X)-\frac{I(A=a)}{e^{*}_{a}(X)}h^{*}_{a}(X)\mid
X\right]\right]+\psi(a)$
$\displaystyle=\operatorname{E}\left[h^{*}_{a}(X)-\frac{1}{e^{*}_{a}(X)}h^{*}_{a}(X)\operatorname{E}\left[I(A=a)\mid
X\right]\right]+\psi(a)$
$\displaystyle=\operatorname{E}\left[h^{*}_{a}(X)-\frac{1}{e^{*}_{a}(X)}h^{*}_{a}(X)\Pr\left[A=a\mid
X\right]\right]+\psi(a)$
$\displaystyle=\operatorname{E}\left[h^{*}_{a}(X)-h^{*}_{a}(X)\right]+\psi(a)$
$\displaystyle=\psi(a).$

Next consider the case when $\widehat{h}_{a}(X)$ is correctly specified, that
is

$h^{*}_{a}(X)=\mathrm{E}\left[L\left(Y,g\left(X^{*}\right)\right)\mid
X,A=a\right]$

and this time we do not make the assumptions that the limit $e^{*}_{a}(X)$ is
equal to $\operatorname{Pr}[A=a\mid X]$. Recall, as shown previously
$\psi(a)=\operatorname{E}\left[\operatorname{E}\left[L(Y,\mu_{\widehat{\beta}}(X^{*}))\mid
X,A=a\right]\right]$.

$\displaystyle\widehat{\psi}_{DR}(a)$ $\displaystyle\stackrel{{\scriptstyle
P}}{{\rightarrow}}\operatorname{E}\left[h^{*}_{a}(X)+\frac{I(A=a)}{e^{*}_{a}(X)}\left(L\left(Y,\mu\left(X^{*}\right)\right)-h^{*}_{a}(X)\right)\right]$
$\displaystyle=\operatorname{E}\left[h^{*}_{a}(X)\right]+\operatorname{E}\left[\frac{I(A=a)}{e^{*}_{a}(X)}\left(L\left(Y,\mu\left(X^{*}\right)\right)-h^{*}_{a}(X)\right)\right]$
$\displaystyle=\psi(a)+\operatorname{E}\left[\frac{I(A=a)}{e^{*}_{a}(X)}\left(L\left(Y,\mu\left(X^{*}\right)\right)-h^{*}_{a}(X)\right)\right]$
$\displaystyle=\psi(a)+\operatorname{E}\left[\operatorname{E}\left[\frac{I(A=a)}{e^{*}_{a}(X)}\left(L\left(Y,\mu\left(X^{*}\right)\right)-h^{*}_{a}(X)\right)\mid
X\right]\right]$
$\displaystyle=\psi(a)+\operatorname{E}\left[\frac{I(A=a)}{e^{*}_{a}(X)}\operatorname{E}\left[\left(L\left(Y,\mu\left(X^{*}\right)\right)-h^{*}_{a}(X)\right)\mid
X\right]\right]$
$\displaystyle=\psi(a)+\operatorname{E}\left[\operatorname{E}\left[\left(L\left(Y,\mu\left(X^{*}\right)\right)-h^{*}_{a}(X)\right)\mid
X,A=a\right]\right]$
$\displaystyle=\psi(a)+\operatorname{E}\left[\operatorname{E}\left[L\left(Y,\mu\left(X^{*}\right)\right)\mid
X,A=a\right]-h^{*}_{a}(X)\right]$
$\displaystyle=\psi(a)+\operatorname{E}\left[h^{*}_{a}(X)-h^{*}_{a}(X)\right]$
$\displaystyle=\psi(a).$

####  C.3.2 Asymptotic distribution

For a random variable $W$ we define notation

$\mathbb{G}_{n}(W)=\sqrt{n}\left(\frac{1}{n}\sum_{i=1}^{n}W_{i}-\mathrm{E}[W]\right).$

and thus the asymptotic representation of $\widehat{\psi}_{DR}(a)$ can be
written

$\displaystyle\sqrt{n}\left(\widehat{\psi}_{DR}(a)-\psi(a)\right)=\mathbb{G}_{n}$
$\displaystyle(H(\widehat{e}_{a}(X),\widehat{h}_{a}(X)))-\mathbb{G}_{n}\left(H\left(e^{*}_{a}(X),h^{*}_{a}(X)\right)\right)$
$\displaystyle+\mathbb{G}_{n}\left(H\left(e^{*}_{a}(X),h^{*}_{a}(X)\right)\right)$
$\displaystyle+\sqrt{n}(\mathrm{E}[H(\widehat{e}_{a}(X),\widehat{h}_{a}(X))]-\psi(a))$

where we add and subtract the term
$\mathbb{G}_{n}\left(H\left(e^{*}_{a}(X),h^{*}_{a}(X)\right)\right)$ and add
another zero term in
$+\sqrt{n}(\mathrm{E}[H(\widehat{e}_{a}(X),\widehat{h}_{a}(X))]-\psi(a))$. For
the first term, Assumption D1 implies

$\mathbb{G}_{n}(H(\widehat{e}_{a}(X),\widehat{h}_{a}(X)))-\mathbb{G}_{n}\left(H\left(e^{*}_{a}(X),h^{*}_{a}(X)\right)\right)=o_{P}(1)$

Let

$Re=\sqrt{n}(\mathrm{E}[H(\widehat{e}_{a}(X),\widehat{h}_{a}(X))]-\psi(a))$

now we have

$\sqrt{n}\left(\widehat{\psi}_{DR}(a)-\psi(a)\right)=\sqrt{n}\left(\frac{1}{n}\sum_{i=1}^{n}\left(H\left(e^{*}_{a}(X_{i}),h^{*}_{a}(X_{i})\right)-\mathrm{E}\left[H\left(e^{*}_{a}(X),h^{*}_{a}(X)\right)\right]\right)\right)+Re+o_{P}(1)$

Let’s try to calculate the upper bound of $Re$. First, note

$n^{-1/2}Re=\underbrace{\mathrm{E}\left[\widehat{h}_{a}(X)\right]}_{R_{1}}+\underbrace{\mathrm{E}\left[\frac{I(A=a)}{\widehat{e}_{a}(X)}\left[L\left(Y,\mu\left(X^{*}\right)\right)-\widehat{h}_{a}(X)\right]\right]}_{R_{2}}-\psi(a).$

We rewrite term $R_{2}$ as:

$\displaystyle R_{2}$
$\displaystyle=\mathrm{E}\left[\frac{I(A=a)}{\widehat{e}_{a}(X)}\left\\{L\left(Y,\mu\left(X^{*}\right)\right)-\widehat{h}_{a}(X)\right\\}\right]$
$\displaystyle=\mathrm{E}\left[\mathrm{E}\left[\frac{I(A=a)}{\widehat{e}_{a}(X)}\left\\{L\left(Y,\mu\left(X^{*}\right)\right)-\widehat{h}_{a}(X)\right\\}\mid
X\right]\right]$
$\displaystyle=\mathrm{E}\left[\frac{1}{\widehat{e}_{a}(X)}\mathrm{E}\left[\frac{I(A=a)}{\operatorname{Pr}[A=a\mid
X]}\operatorname{Pr}[A=a\mid
X]\left\\{L\left(Y,\mu\left(X^{*}\right)\right)-\widehat{h}_{a}(X)\right\\}\mid
X\right]\right]$
$\displaystyle=\mathrm{E}\left[\frac{1}{\widehat{e}_{a}(X)}\mathrm{E}\left[\operatorname{Pr}[A=a\mid
X]\left\\{L\left(Y,\mu\left(X^{*}\right)\right)-\widehat{h}_{a}(X)\right\\}\mid
X,A=a\right]\right]$
$\displaystyle=\mathrm{E}\left[\frac{1}{\widehat{e}_{a}(X)}\operatorname{Pr}[A=a\mid
X]\left\\{\mathrm{E}\left[L\left(Y,\mu\left(X^{*}\right)\right)\mid
X,A=a\right]-\widehat{h}_{a}(X)\right\\}\right]$

Combining the above gives

$\displaystyle n^{-1/2}Re$
$\displaystyle=\mathrm{E}\left[\widehat{h}_{a}(X)\right]+\mathrm{E}\left[\frac{I(A=a)}{e_{a}^{\prime}(X)}\left[L\left(Y,\mu\left(X^{*}\right)\right)-h_{a}^{\prime}(X)\right]\right]-\psi(a)$
$\displaystyle=\mathrm{E}\left[\widehat{h}_{a}(X)\right]+\mathrm{E}\left[\frac{1}{\widehat{e}_{a}(X)}\operatorname{Pr}[A=a\mid
X]\left\\{\mathrm{E}\left[L\left(Y,\mu\left(X^{*}\right)\right)\mid
X,A=a\right]-\widehat{h}_{a}(X)\right\\}\right]$
$\displaystyle-\operatorname{E}\left[\operatorname{E}\left[L(Y,\mu_{\widehat{\beta}}(X^{*}))\mid
X,A=a\right]\right]$
$\displaystyle=\operatorname{E}\left[\left\\{\mathrm{E}\left[L\left(Y,\mu\left(X^{*}\right)\right)\mid
X,A=a\right]-\widehat{h}_{a}(X)\right\\}\times\left\\{\frac{1}{\widehat{e}_{a}(X)}\operatorname{Pr}[A=a\mid
X]-1\right\\}\right]$

Using the Cauchy-Schwartz inequality we get.

$\displaystyle Re$
$\displaystyle\leq\sqrt{n}\left(\mathrm{E}\left[\left\\{\mathrm{E}\left[L\left(Y,\mu\left(X^{*}\right)\right)\mid
X,A=a\right]-\widehat{h}_{a}(X)\right\\}^{2}\right]\right)^{1/2}$
$\displaystyle\times\left(\mathrm{E}\left[\left\\{\frac{1}{\widehat{e}_{a}(X)}\operatorname{Pr}[A=a\mid
X]-1\right\\}^{2}\right]\right)^{1/2}$ $\displaystyle\leq$
$\displaystyle\sqrt{n}O_{P}\left(\left\|\mathrm{E}\left[L\left(Y,\mu\left(X^{*}\right)\right)\mid
X,A=a\right]-\widehat{h}_{a}(X)\right\|_{2}^{2}\times\Big{\|}\widehat{e}_{a}(X)-\operatorname{Pr}[A=a\mid
X]\Big{\|}_{2}^{2}\right)$

If both models $\widehat{e}_{a}(X)$ and $\widehat{h}_{a}(X)$ are correctly
specified and converge at a combined rate faster than $\sqrt{n}$, then
$Re=o_{P}(1)$ and

$\displaystyle\sqrt{n}\left(\widehat{\psi}_{DR}(a)-\psi(a)\right)$
$\displaystyle=\sqrt{n}\left(\frac{1}{n}\sum_{i=1}^{n}H\left(\operatorname{Pr}\left[A=a\mid
X_{i}\right],\mathrm{E}\left[L\left(Y,g\left(X^{*}\right)\right)\mid
A=a,X_{i}\right]\right)\right.$
$\displaystyle\left.-\mathrm{E}\left[H\left(\operatorname{Pr}[A=a\mid
X],\mathrm{E}\left[L\left(Y,g\left(X^{*}\right)\right)\mid
A=a,X\right]\right)\right]\right)+o_{P}(1)$

By the central limit theorem,

$\sqrt{n}\left(\frac{1}{n}\sum_{i=1}^{n}H\left(e^{*}_{a}(X_{i}),h^{*}_{a}(X_{i})\right)-\mathrm{E}\left[H\left(e^{*}_{a}(X_{i}),h^{*}_{a}(X_{i})\right)\right]\right)\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N\left(0,\operatorname{Var}\left[H\left(e^{*}_{a}(X),h^{*}_{a}(X)\right)\right]\right)$

completing the proof.

##  Appendix D Risk calibration curve

When the outcome is binary, another common metric of model performance is risk
calibration. Calibration is a measure of the relibability of the risk
estimates produced by the fitted model $\mu_{\widehat{\beta}}(X^{*})$. For
instance, among a sample of patients who receive a risk prediction of 17% does
the outcome really occur for roughly 17% of them over the follow up period?
This can be nonparametrically evalutated across a range of risks by estimating
the so-called “calibration” curve, i.e. the observed risk as a function of the
predicted risk. For counterfactual predictions the relevant calibration curve
is the counterfactual risk that would be observed under intervetion $A=a$ as a
function of the predicted risk, or

$\psi_{\widehat{\beta}}(a)=\operatorname{E}[I(Y^{a}=1)\mid\mu_{\widehat{\beta}}(x^{*})].$
(A9)

###  D.1 Identification

Here we show that the counterfactual calibration curve
$\psi_{\widehat{\beta}}(a)=\operatorname{E}[I(Y^{a}=1)\mid\mu_{\widehat{\beta}}(x^{*})]$
is identified using the observed data under the assumptions of section 4.

###### Theorem 5

Under conditions A1-A3, the risk calibration curve is identified by the
observed data functionals

$\psi_{\widehat{\beta}}(a)=\operatorname{E}[\operatorname{E}\\{I(Y=1)\mid
X,A=a,\mu_{\widehat{\beta}}(x^{*}),D_{test}=1\\}\mid\mu_{\widehat{\beta}}(x^{*}),D_{test}=1]$
(A10)

and

$\psi_{\widehat{\beta}}(a)=\operatorname{E}\left[\frac{I(A=a)}{\Pr[A=a\mid
X,\mu_{\widehat{\beta}}(x^{*}),D_{test}=1]}I(Y=1)\mid\mu_{\widehat{\beta}}(x^{*}),D_{test}=1\right]$
(A11)

in the test set.

Proof. For the first representation we have

$\displaystyle\psi_{\widehat{\beta}}(a)$
$\displaystyle=\operatorname{E}[I(Y^{a}=1)\mid\mu_{\widehat{\beta}}(x^{*})]$
$\displaystyle=\operatorname{E}[I(Y^{a}=1)\mid\mu_{\widehat{\beta}}(x^{*}),D_{test}=1]$
$\displaystyle=\operatorname{E}[\operatorname{E}\\{I(Y^{a}=1)\mid
X,\mu_{\widehat{\beta}}(x^{*}),D_{test}=1\\}\mid\mu_{\widehat{\beta}}(x^{*}),D_{test}=1]$
$\displaystyle=\operatorname{E}[\operatorname{E}\\{I(Y^{a}=1)\mid
X,A=a,\mu_{\widehat{\beta}}(x^{*}),D_{test}=1\\}\mid\mu_{\widehat{\beta}}(x^{*}),D_{test}=1]$
$\displaystyle=\operatorname{E}[\operatorname{E}\\{I(Y=1)\mid
X,A=a,\mu_{\widehat{\beta}}(x^{*}),D_{test}=1\\}\mid\mu_{\widehat{\beta}}(x^{*}),D_{test}=1]$

where the first line follows from the definition of
$\psi_{\widehat{\beta}}(a)$, the second from random sampling of the test set,
the third from the law of iterated expectations, the fourth from the
exchangeability condition A1, and the fifth from the consistency condition A2.
Recall that $X^{*}$ is a subset of $X$. For the second representation, we show
that it is equivalent to the first

$\displaystyle\psi_{\widehat{\beta}}(a)$
$\displaystyle=\operatorname{E}[\operatorname{E}\\{I(Y=1)\mid
X,A=a,\mu_{\widehat{\beta}}(x^{*}),D_{test}=1\\}\mid\mu_{\widehat{\beta}}(x^{*}),D_{test}=1]$
$\displaystyle=\operatorname{E}\left[\operatorname{E}\left\\{\frac{I(A=a)}{\Pr[A=a\mid
X,\mu_{\widehat{\beta}}(x^{*}),D_{test}=1]}I(Y=1)\mid
X,\mu_{\widehat{\beta}}(x^{*}),D_{test}=1\right\\}\mid\mu_{\widehat{\beta}}(x^{*}),D_{test}=1\right]$
$\displaystyle=\operatorname{E}\left[\frac{I(A=a)}{\Pr[A=a\mid
X,\mu_{\widehat{\beta}}(x^{*}),D_{test}=1]}I(Y=1)\mid\mu_{\widehat{\beta}}(x^{*}),D_{test}=1\right]$

where the second line follows from the definition of conditional expectation,
and the last reverses the law of iterated expectations.

###  D.2 Estimation

Unlike previous sections, estimation of the full risk calibration curve using
sample analogs of the identified expressions A10 and A11 is generally
infeasible because they are conditional on a continuous risk score. Instead
analysts typically perform either kernel or binned estimation of the
calibration curve functional. In the case of the counterfactual risk
calibration curve under a hypothetical intervention, the expression above
suggests modifying these approaches either through the use of inverse
probability weights or an outcome model.

##  Appendix E Area under ROC curve

Another common metric for the performance of a risk prediction model
$\mu_{\beta}(X^{*})$ is the area under the receiver operating characteristic
(ROC) curve, often referred to as simply the area under the curve (AUC). The
AUC can be interpreted as the probability that a randomly sampled observation
with the outcome has a higher predicted value than a randomly sampled
observation without the outcome. In that sense, it is a measure of the
discriminative ability of the model, i.e. the ability to distinguish between
cases and noncases. For counterfactual predictions the relevant AUC is the
counterfactual AUC that would be observed under intervetion $A=a$, or

$\psi_{\widehat{\beta}}(a)=\operatorname{E}[I\left(\mu_{\beta}(X^{*}_{i})>\mu_{\beta}(X^{*}_{j})\right)\mid
Y_{i}^{a}=1,Y_{j}^{a}=0].$ (A12)

###  E.1 Identification

Here, we show that the counterfactual AUC $\psi_{\widehat{\beta}}(a)$ is
identified by the observed data under a modified set of identification
conditions, namely:

  1. E1.

Exchangeability. $Y^{a}\perp\\!\\!\\!\perp A\mid X$

  2. E2.

Consistency. $Y^{a}=Y$ if $A=a$

  3. E3.

Positivity. (i) $\Pr(A=a|X=x)>0$ for all $x$ that have positive density in
$f(X,A=a)$, (ii)
$\mathrm{E}\left[\Pr[Y=1|X_{i},A=a]\Pr[Y=0|X_{j},A=a]\right]>0$, where $i$ is
a random observation that has the outcome and $j$ is random observation
without the outcome.

###### Theorem 6

Under conditions E1-E3, the counterfactual AUC is identified by the observed
data functionals in the test set

$\psi_{\widehat{\beta}}(a)=\frac{\mathrm{E}\left[I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*})\right)m_{a}(X_{i},X_{j})\right]}{\mathrm{E}\left[m_{a}(X_{i},X_{j})\right]}$
(A13)

and

$\psi_{\widehat{\beta}}(a)=\frac{\mathrm{E}\left[\frac{I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*}),Y_{i}=1,Y_{j}=0,A_{i}=a,A_{j}=a\right)}{\pi_{a}(X_{i},X_{j})}\right]}{\mathrm{E}\left[\frac{I\left(Y_{i}=1,Y_{j}=0,A_{i}=a,A_{j}=a\right)}{\pi_{a}(X_{i},X_{j})}\right]}$
(A14)

where the subscripts $i$ and $j$ denote a random pair of observations from the
test set where

$m_{a}(X_{i},X_{j})=\operatorname{Pr}\left[Y_{i}=1\mid
X_{i},A_{i}=a,D_{test,i}=1\right]\Pr\left[Y_{j}=0\mid
X_{j},A_{j}=a,D_{test,j}=1\right]$

and

$\pi_{a}(X_{i},X_{j})=\Pr\left[A_{i}=a\mid
X_{i},D_{test,i}=1\right]\Pr\left[A_{j}=a\mid X_{j},D_{test,j}=1\right]$

for a pair of covariate vectors $X_{i}$ and $X_{j}$.

Proof. For the first representation we have

$\displaystyle\psi_{\widehat{\beta}}(a)$
$\displaystyle=\mathrm{E}\left[I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*})\right)\mid
Y^{a}_{i}=1,Y^{a}_{j}=0\right]$
$\displaystyle=\frac{\mathrm{E}\left[I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*}),Y^{a}_{i}=1,Y^{a}_{j}=0\right)\right]}{\operatorname{Pr}\left[Y^{a}_{i}=1,Y^{a}_{j}=0\right]}$
$\displaystyle=\frac{\mathrm{E}\left[\mathrm{E}\left[I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*}),Y^{a}_{i}=1,Y^{a}_{j}=0\right)\mid
X_{i},X_{j}\right]\right]}{\mathrm{E}\left[\operatorname{Pr}\left[Y^{a}_{i}=1,Y^{a}_{j}=0\mid
X_{i},X_{j}\right]\right]}$
$\displaystyle=\frac{\mathrm{E}\left[I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*})\right)\operatorname{Pr}\left[Y^{a}_{i}=1,Y^{a}_{j}=0\mid
X_{i},X_{j}\right]\right]}{\mathrm{E}\left[\operatorname{Pr}\left[Y^{a}_{i}=1,Y^{a}_{j}=0\mid
X_{i},X_{j}\right]\right]}$
$\displaystyle=\frac{\mathrm{E}\left[I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*})\right)\operatorname{Pr}\left[Y^{a}_{i}=1,Y^{a}_{j}=0\mid
X_{i},X_{j},A_{i}=a,A_{j}=a\right]\right]}{\mathrm{E}\left[\operatorname{Pr}\left[Y^{a}_{i}=1,Y^{a}_{j}=0\mid
A_{i}=a,A_{j}=a,X_{i},X_{j}\right]\right]}$
$\displaystyle=\frac{\mathrm{E}\left[I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*})\right)\operatorname{Pr}\left[Y^{a}_{i}=1\mid
X_{i},A_{i}=a\right]\Pr\left[Y^{a}_{j}=0\mid
X_{j},A_{j}=a\right]\right]}{\mathrm{E}\left[\operatorname{Pr}\left[Y^{a}_{i}=1\mid
X_{i},A_{i}=a\right]\Pr\left[Y^{a}_{j}=0\mid X_{j},A_{j}=a\right]\right]}$
$\displaystyle=\frac{\mathrm{E}\left[I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*})\right)\operatorname{Pr}\left[Y_{i}=1\mid
X_{i},A_{i}=a\right]\Pr\left[Y_{j}=0\mid
X_{j},A_{j}=a\right]\right]}{\mathrm{E}\left[\operatorname{Pr}\left[Y_{i}=1\mid
X_{i},A_{i}=a\right]\Pr\left[Y_{j}=0\mid X_{j},A_{j}=a\right]\right]}$
$\displaystyle=\frac{\mathrm{E}\left[I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*})\right)\operatorname{Pr}\left[Y_{i}=1\mid
X_{i},A_{i}=a\right]\Pr\left[Y_{j}=0\mid
X_{j},A_{j}=a\right]\right]}{\mathrm{E}\left[\operatorname{Pr}\left[Y_{i}=1\mid
X_{i},A_{i}=a\right]\Pr\left[Y_{j}=0\mid X_{j},A_{j}=a\right]\right]}$
$\displaystyle=\frac{\mathrm{E}\left[I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*})\right)m_{a}(X_{i},X_{j})\right]}{\mathrm{E}\left[m_{a}(X_{i},X_{j})\right]}$

where the first line follows from the definition of
$\psi_{\widehat{\beta}}(a)$, the second from the definition of conditional
probability, the third from the law of iterated expectations, the fourth from
the definition of conditional expectation, the fifth from the exchangeability
condition E1, the sixth from independence of potential outcomes, the seventh
from the consistency condition E2, the eighth from random sampling of the test
set, and the ninth applies the definition of $m_{a}(X_{i},X_{j})$. Recall that
$X^{*}$ is a subset of $X$. For the second representation, we will show that
it is equivalent to the first. Starting from line five above

$\displaystyle\psi_{\widehat{\beta}}(a)$
$\displaystyle=\frac{\mathrm{E}\left[I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*})\right)\operatorname{Pr}\left[Y^{a}_{i}=1,Y^{a}_{j}=0\mid
X_{i},X_{j},A_{i}=a,A_{j}=a\right]\right]}{\mathrm{E}\left[\operatorname{Pr}\left[Y^{a}_{i}=1,Y^{a}_{j}=0\mid
A_{i}=a,A_{j}=a,X_{i},X_{j}\right]\right]}$
$\displaystyle=\frac{\mathrm{E}\left[I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*})\right)\operatorname{Pr}\left[Y_{i}=1,Y_{j}=0\mid
X_{i},X_{j},A_{i}=a,A_{j}=a\right]\right]}{\mathrm{E}\left[\operatorname{Pr}\left[Y_{i}=1,Y_{j}=0\mid
A_{i}=a,A_{j}=a,X_{i},X_{j}\right]\right]}$
$\displaystyle=\frac{\mathrm{E}\left[I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*})\right)\frac{\operatorname{Pr}\left[Y_{i}=1,Y_{j}=0,A_{i}=a,A_{j}=a\mid
X_{i},X_{j}\right]}{\operatorname{Pr}\left[A_{i}=a,A_{j}=a\mid
X_{i},X_{j}\right]}\right]}{\mathrm{E}\left[\frac{\operatorname{Pr}\left[Y_{i}=1,Y_{j}=0,A_{i}=a,A_{j}=a\mid
X_{i},X_{j}\right]}{\operatorname{Pr}\left[A_{i}=a,A_{j}=a\mid
X_{i},X_{j}\right]}\right]}$
$\displaystyle=\frac{\mathrm{E}\left[\mathrm{E}\left[I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*})\right)\frac{\operatorname{Pr}\left[Y_{i}=1,Y_{j}=0,A_{i}=a,A_{j}=a\mid
X_{i},X_{j}\right]}{\operatorname{Pr}\left[A_{i}=a,A_{j}=a\mid
X_{i},X_{j}\right]}\mid
X_{i},X_{j}\right]\right]}{\mathrm{E}\left[\mathrm{E}\left[\frac{\operatorname{Pr}\left[Y_{i}=1,Y_{j}=0,A_{i}=a,A_{j}=a\mid
X_{i},X_{j}\right]}{\operatorname{Pr}\left[A_{i}=a,A_{j}=a\mid
X_{i},X_{j}\right]}\mid X_{i},X_{j}\right]\right]}$
$\displaystyle=\frac{\mathrm{E}\left[\frac{I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*})\right)}{\operatorname{Pr}\left[A_{i}=a\mid
X_{i}\right]\operatorname{Pr}\left[A_{j}=a\mid
X_{j}\right]}\operatorname{Pr}\left[Y_{i}=1,Y_{j}=0,A_{i}=a,A_{j}=a\mid
X_{i},X_{j}\right]\right]}{\mathrm{E}\left[\frac{\operatorname{Pr}\left[Y_{i}=1,Y_{j}=0,A_{i}=a,A_{j}=a\mid
X_{i},X_{j}\right]}{\operatorname{Pr}[A_{i}=a\mid X_{i}]\Pr[A_{j}=a\mid
X_{j}]}\right]}$
$\displaystyle=\frac{\mathrm{E}\left[\frac{I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*}),Y_{i}=1,Y_{j}=0,A_{i}=a,A_{j}=a\right)}{\operatorname{Pr}\left[A_{i}=a\mid
X_{i}\right]\operatorname{Pr}\left[A_{j}=a\mid
X_{j}\right]}\right]}{\mathrm{E}\left[\frac{I\left(Y_{i}=1,Y_{j}=0,A_{i}=a,A_{j}=a\right)}{\operatorname{Pr}[A_{i}=a\mid
X_{i}]\Pr[A_{j}=a\mid X_{j}]}\right]}$
$\displaystyle=\frac{\mathrm{E}\left[\frac{I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*}),Y_{i}=1,Y_{j}=0,A_{i}=a,A_{j}=a\right)}{\pi_{a}(X_{i},X_{j})}\right]}{\mathrm{E}\left[\frac{I\left(Y_{i}=1,Y_{j}=0,A_{i}=a,A_{j}=a\right)}{\pi_{a}(X_{i},X_{j})}\right]}$

where the second line follows from consistency $E2$, the third from the
definition of conditional probability, the fourth from iterated expectations,
the fifth removes the constant fraction outside expectation, the sixth
reverses the law of iterated expectations and the last applies random sampling
of the test set and the definition of $\pi_{a}(X_{i},X_{j})$.

###  E.2 Plug-in estimation

Using sample analogs for the identified expressions A13 and A14, we obtain two
plug-in estimators for the counterfactual AUC

$\widehat{\psi}_{OM}(a)=\frac{\sum_{i\neq
j}^{n}\widehat{h}_{a}(X_{i})(1-\widehat{h}_{a}(X_{j}))I(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*}),D_{test,i}=1,D_{test,j}=1)}{\sum_{i\neq
j}^{n}\widehat{h}_{a}(X_{i})(1-\widehat{h}_{a}(X_{j}))I(D_{test,i}=1,D_{test,j}=1)}$

and

$\widehat{\psi}_{IPW}(a)=\frac{{{\sum}}_{i\neq
j}^{n}\dfrac{I\left(\mu_{\widehat{\beta}}(X_{i}^{*})>\mu_{\widehat{\beta}}(X_{j}^{*}),Y_{i}=1,Y_{j}=0,A_{i}=a,A_{j}=a,D_{test,i}=1,D_{test,j}=1\right)}{\widehat{e}_{a}(X_{i})\widehat{e}_{a}(X_{j})}}{{{\sum}}_{i\neq
j}^{n}\dfrac{I\left(Y_{i}=1,Y_{j}=0,A_{i}=a,A_{j}=a,D_{test,i}=1,D_{test,j}=1\right)}{\widehat{e}_{a}(X_{i})\widehat{e}_{a}(X_{j})}}$

where
$\widehat{m}_{a}(X_{i},X_{j})=\widehat{h}_{a}(X_{i})(1-\widehat{h}_{a}(X_{j}))$
and $\widehat{h}_{a}$ is an estimator for
$\operatorname{Pr}[Y_{i}=1|X_{i},A_{i}=a,D_{test,i}=1]$ and where
$\widehat{\pi}_{a}(X_{i},X_{j})=\widehat{e}_{a}(X_{i})\widehat{e}_{a}(X_{j})$
and $\widehat{e}_{a}$ is an estimator for $\Pr[A_{i}=a|X_{i},D_{test,i}=1]$.
Here, we call the first plug-in estimator the outcome model estimator
$\widehat{\psi}_{OM}$ and the second the inverse probability weighted
estimator $\widehat{\psi}_{IPW}$.

##  Appendix F Additional application details

The Multi-Ethnic Study on Atherosclerosis (MESA) study is a population-based
sample of 6,814 men and women aged 45 to 84 drawn from six communities
(Baltimore; Chicago; Forsyth County, North Carolina; Los Angeles; New York;
and St. Paul, Minnesota) in the United States between 2000 and 2002. The
sampling procedure, design, and methods of the study have been described
previously [31]. Study teams conducted five examination visits between 2000
and 2011 in 18 to 24 month intervals focused on the prevalence, correlates,
and progression of subclinical cardiovascular disease. These examinations
included assessments of lipid-lowering (primarily statins) and other
medication use as well as cardiovascular risk factors such as systolic blood
pressure, serum cholesterol, cigarette smoking, height, weight, and diabetes.

Our goal was to emulate a single-arm trial corresponding to the AHA guidelines
on initiation of statin therapy for primary prevention of cardiovascular
disease in the MESA cohort and use the emulated trial to develop a prediction
model for the treatment-naive risk. The AHA guidelines stipulate that patients
aged 40 to 75 with serum LDL cholesterol levels between 70 mg/dL and 190 mg/dL
and no history of cardiovascular disease should initiate statins if their risk
exceeds 7.5%. Therefore, we considered MESA participants who completed the
baseline examination, had no recent history of statin use, no history of
cardiovascular disease, and who met the criteria described in the guidelines
(excluding the risk threshold) as eligible to participate in the trial. The
primary endpoint was time to atherosclerotic cardiovascular disease (ASCVD),
defined as nonfatal myocardial infarction, coronary heart disease death, or
ischemic stroke.

Follow up began at the second examination cycle to enable a “wash out” period
for statin use and to ensure adequate pre-treatment covariates to control
confouding. We constructed a sequence of nested trials starting at each
examination cycle from exam 2 through exam 5 and pooled the results from all 4
trials into a single analysis and used a robust variance estimator to account
for correlation among duplicated participants. In each nested trial, we used
the corresponding questionnaire to determine eligibility as well as statin
initiators versus non-initiators. Because the exact timing of statin
initiation was not known with precision, in each trial, we estimated the start
of follow up for initiators and non-initators by drawing a random month
between their current and previous examinations. We explored alternative
definitions of the start of follow up in sensitivity analyses in the appendix.
To mimic the targeted single-arm trial we limited to non-initiators for
development of the prediction models.

###  F.1 Propensity score models

In the emulated single arm trial, statin initiation can be viewed as “non-
adherence” which can be adjusted for by inverse probability weighting,
therefore we censored participants when they initiated statins. To calculate
the weights, we estimated two logistic regression models: one for the
probability of remaining untreated given past covariate history (denominator
model) and one for probability of remaining untreated given the selected
baseline predictors (numerator model). In the denominator model we included
the following covariates:

  * •

Demographic factors \- Age, gender, marital status, education, race/ethnicity,
employment, health insurance status, depression, perceived discrimination,
emotional support, anger and anxiety scales, and neighborhood score.

  * •

Risk factors \- Systolic and diastolic blood pressure, serum cholesterol
levels (LDL, HDL, Triglycerides), hypertension, diabetes, waist circumference,
smoking, alcohol consumption, exercise, family history of CVD, calcium score,
hypertrophy on ECG, CRP, IL-6, number of pregnancies, oral contraceptive use,
age of menopause.

  * •

Medication use \- Anti-hypertensive use, insulin use, daily aspirin use, anti-
depressant use, vasodilator use, anti-arryhtmic use.

Time-varying demographic factors and risk factors were lagged such that values
from the previous examination cycle were used.