###### Proof.

.: 5pt

######  Robust Estimation of Loss-Based Measures of Model Performance under
Covariate Shift

###### Abstract

We present methods for estimating loss-based measures of the performance of a
prediction model in a target population that differs from the source
population in which the model was developed, in settings where outcome and
covariate data are available from the source population but only covariate
data are available on a simple random sample from the target population. Prior
work adjusting for differences between the two populations has used various
weighting estimators with inverse odds or density ratio weights. Here, we
develop more robust estimators for the target population risk (expected loss)
that can be used with data-adaptive (e.g., machine learning-based) estimation
of nuisance paramaters. We examine the large-sample properties of the
estimators and evaluate finite sample performance in simulations. Last, we
apply the methods to data from lung cancer screening using nationally
representative data from the National Health and Nutrition Examination Survey
(NHANES) and extend our methods to account for the complex survey design of
the NHANES.

Keywords: transportability, covariate shift, domain adaptation, MSE, double
robustness, weighting

##  1 Introduction

Ideally, a prediction model should be evaluated using data from the target
population where it will be applied, but typically the source data used for
model building are not obtained from a random sample of that target population
[1] and model performance in the source data may not reflect performance in
the target population. A major reason why model performance estimates may not
transport to the target population is “covariate shift,” that is, the presence
of differences in the covariate distribution between the population underlying
the source data (i.e., the source population) and the target population [2, 3,
4]. Under covariate shift, the conditional distribution of the outcome given
the covariates is the same in both the source and target population, but the
covariate distributions of the two populations are different (i.e., the
populations have a different “case-mix”). Such differences affect model
performance as evaluated by measures that are averages of a loss function over
the target population distribution (e.g., the mean squared error, MSE; mean
absolute error; the Brier score [5]).

To understand how covariate shift affects model performance, it is useful to
consider “prediction error modifiers” – covariates that are associated with
prediction error, as assessed by a specific loss function, for a given model
[6]. When prediction error modifiers have a different distribution between the
source and target population, estimators of loss-based measures of model
performance that only use source data are biased for target population model
performance [7]. Previous efforts to correct this bias used importance-
weighting approaches [8] to re-weight observations in the source data by the
ratio of the covariate density of the target and source population or,
equivalently, by the inverse of the conditional odds of being from the source
population [6], in order to construct asymptotically unbiased estimators of
the target population risk (i.e., the expected loss in the target population).
The conditional odds (or the covariate density ratio) are almost always
unknown and need to be estimated using statistical models. For the importance-
weighting estimators to consistently estimate the target population risk, the
models for the conditional odds needs to be correctly specified.

Even when the model for the conditional odds of being from the source
population is correctly specified, importance-weighting estimators are
inefficient [9]. Furthermore, asymptotically valid inference for importance-
weighting estimators requires estimators of the conditional odds to converge
at $\sqrt{n}$ rate [10]. That precludes using data-adaptive estimators (e.g.,
machine learning estimators) of the density ratio or the conditional odds
because such estimators converge at slower than $\sqrt{n}$ rate. Data-adaptive
estimators, however, are very appealing in applied work because subject matter
knowledge is typically inadequate to determine the correct specification of
the conditional odds models, particularly when the covariates that differ in
distribution between the source and target population are high-dimensional or
have multiple continuous components.

In this paper, we develop a doubly robust estimator for the target population
risk that involves estimating both the expected loss conditional on covariates
and the probability of participation in the source. Our estimator is
consistent for the target population risk if at least one of these models is
correctly specified, but not necessarily both, and can be used for
asymptotically valid inference even if the models are estimated using methods
that converge at a rate slower than $\sqrt{n}$ (i.e., allowing the use of
data-adaptive estimation). In the process of developing the doubly robust
estimator, we also develop a novel conditional loss modeling-based estimator
that relies on estimating the expected loss conditional on covariates. We
provide identifiability conditions, identification results, and large-sample
properties for the doubly robust estimator. We compare the finite-sample
performance of the doubly robust and conditional loss modeling-based
estimators against importance-weighting estimators in simulations. Last, we
apply the methods to estimate model performance of a prediction model for lung
cancer diagnosis built using data from the National Lung Screening Trial
(NLST) in a target population of people eligible for lung cancer screening in
the US. The sample from the target population comes from the National Health
and Nutrition Examination Survey (NHANES) a complex survey that involves
multi-stage clustering and variable probability sampling. We show how to
modify our estimators to account for the complex sampling design and to
incorporate weights that account for the oversampling of certain subgroups,
survey non-response, and post-stratification adjustments.

##  2 Identificaton of the target population risk

###  2.1 Setup, study design, and targets of inference

Let $Y$ be an outcome and $X$ a covariate vector. Let $X^{*}$ be a vector that
contains a subset of the covariates in $X$, and $g(X^{*})$ a prediction model
for the conditional expectation of $Y$ given $X^{*}$ that we are interested in
evaluating the performance of in a target population of substantive interest.
We consider a setting where we have access to outcome and covariate data on a
sample from the source population $\\{(X_{i},Y_{i}):i=1,\ldots,n_{1}\\}$ and
covariate data on a separately obtained random sample from the target
population $\\{X_{i}:i=1,\ldots,n_{0}\\}$. Let $D$ be an indicator of being
from the source population (i.e., $D=1$ if an observation comes from the
source population and $D=0$ if an observation comes from the target
population). In this setup, the data used to estimate model performance in the
target population is the combined source population sample and the target
population sample $\\{(X_{i},D_{i},D_{i}\times
Y_{i}):i=1,\ldots,n=n_{1}+n_{0}\\}$. We also assume that the model is fit
using data that is independent of the data used to assess model performance
(e.g., the model is an external model or it is fit using training data from
the source population and an independent set of test data from the source
population is used for model assessment).

We focus on etimation of loss-based measures of model performance because many
popular evaluation measures are loss-based (e.g., mean squared error, Brier
loss, and the absolute loss). A loss function $L(Y,g(X^{*}))$ measures the
discrepancy between the observed outcome $Y$ and a prediction of the outcome
$g(X^{*})$. We define the target parameter of interest as the target
population risk (i.e., expected loss in the target population):

$\psi=\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|D=0].$

Importantly, we _do not_ assume that the model $g(X^{*})$ is correctly
specified when developing our estimator, which implies that we can estimate
model performance in the target population for both correctly specified and
misspecified models as long as the identifiability assumptions listed below
hold.

###  2.2 Identification analysis

#### Identifiability conditions:

We will argue that $\psi$ is identifiable, that is, can be written as a
functional of the observed data distribution, under the following conditions
which are fairly standard in the literature on transporting prediction models
to a target population [7, 8]:

  * A1.

Conditional exchangeability over $D$: For all $x$ such that $f(X=x,D=0)\neq
0$,

$F(Y|X=x,D=0)=F(Y|X=x,D=1),$

where $F(Y|X=x,D=0)$ and $F(Y|X=x,D=1)$ are the cumulative distribution
functions in the target and source population, respectively, and $f(X,D=0)$ is
the covariate density in the target population.

  * A2.

Positivity of participation in the source: For all $x$ such that
$f(X=x,D=0)\neq 0$,

$\Pr[D=1|X=x]>0.$

Assumption A1 implies that the joint distribution of $(X,Y)$ in the source and
target population differ only in terms of the covariate distribution (i.e.,
the setting of covariate shift); however, we _do not_ assume that the
conditional distribution of $Y$ given $X$ is known or can be consistently
estimated. Assumption A2 informally says that all covariate patterns in the
target population have a positive probability of appearing in the source data.

#### Identification:

We now give an identification result for the target population model
performance, under conditions A1 and A2 [6].

######  Theorem 1.

If conditions A1 and A2 hold, then $\psi$ can be written as the observed data
functional

$\psi=\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|D=0]=\operatorname{\textnormal{\mbox{E}}}[\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]|D=0];$
(1)

or, using a weighting representation,

$\psi=\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}\left[\frac{I(D=1)\Pr[D=0|X]}{\Pr[D=1|X]}L(Y,g(X^{*}))\right].$
(2)

We assume that data from the source and target population are separately
sampled with unknown but typically expected to be unequal sampling fractions
in a “non-nested design”. In this situation the probabilities $\Pr[D=0|X]$ and
$\Pr[D=0]$ appearing in expression (2) are not identifiable (without
additional information) [11, 6] since under the non-nested design the ratio
$n_{0}/(n_{0}+n_{1})$ does not necessarily converge to the population
probability of being in the target population, that is $\Pr[D=0]$.
Nevertheless, as we discuss next, the target parameter $\psi$ is still
identifiable under this biased sampling design.

#### Identification under the biased sampling design:

Let $R$ be an indicator for sampling under the non-nested design (i.e., $R=1$
if an observation is sampled from the source or target population; otherwise
$R=0$). Because the data available are separately obtained random samples from
the source population ($D=1$) and from the target population ($D=0$), we have
the following theorem [6]:

######  Theorem 2.

Under the biased sampling design, the risk in the target population is
identifiable by the observed data functional

$\psi=\operatorname{\textnormal{\mbox{E}}}[\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1,R=1]|D=0,R=1].$
(3)

Furthermore, the right-hand-side of the above equation can be re-expressed as

$\psi=\frac{1}{\Pr[D=0|R=1]}\operatorname{\textnormal{\mbox{E}}}\left[\frac{I(D=1)\Pr[D=0|X,R=1]}{\Pr[D=1|X,R=1]}L(Y,g(X^{*}))\Bigg{|}R=1\right].$
(4)

Because all quantities in expressions (3) and (4) condition on the observed
data ($R=1$), the target parameter $\psi$ is identifiable under the biased
sampling design and we can interpret the expectations in Theorem 1 as being
taken with respect to densities from the biased sampling design. An
alternative design is a nested design where the source population is nested
within a larger target population [11, 12]; the results presented here can be
extended to nested designs.

##  3 Estimation of the target population risk

###  3.1 Sample-analog estimators

Using sample analogs in expression (1) we obtain the following conditional
loss-based estimator:

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

where $\widehat{h}(X)$ is an estimator for
$\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$. For simplicity of
notation we suppress the dependency of $\widehat{h}(X)$ on $g(X^{*})$ and the
loss function $L(\cdot,\cdot)$. The conditional loss estimator
$\widehat{\psi}_{CL}$ is a consistent estimator for $\psi$ if the model
$\widehat{h}(X)$ is correctly specified (i.e., $\widehat{h}(X)$ converges in
probability to $\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$).

Using sample analogs in expression (2) we obtain the following weighting
estimator [13, 14, 6]:

$\widehat{\psi}_{IW}=\frac{1}{n_{0}}\sum_{i=1}^{n}\frac{I(D_{i}=1)(1-\widehat{p}(X_{i}))}{\widehat{p}(X_{i})}L(Y_{i},g(X^{*}_{i})),$
(6)

where $\widehat{p}(X)$ is an estimator for the conditional probability of
belonging to the source population, $\Pr[D=1|X]$. The weighting estimator
$\widehat{\psi}_{IW}$ is a consistent estimator for $\widehat{\psi}$ if the
model $\widehat{p}(X)$ is correctly specified (i.e., $\widehat{p}(X)$
converges in probability to $\Pr[D=1|X]$). An advantage of
$\widehat{\psi}_{IW}$ compared to the importance weighting estimator proposed
by [8] is that it depends on the ratio of conditional probabilities of being
in the source and the target population rather than the ratio of the source
and target population densities. There are many established methods for
estimating conditional probabilities with high dimensional covariates [15],
whereas density estimation becomes challenging even with moderately
dimensional covariate vectors.

###  3.2 Doubly robust estimator

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 being in the source data. Both
estimators also inherit the rate of convergence of the estimator used to
estimate the nuisance parameter needed for their implementation ($\Pr[D=1|X]$
or $\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$). Hence if the
conditional loss or weighting estimators are combined with data-adaptive
estimators for $\widehat{p}(X)$ or $\widehat{h}(X)$ that converge at a rate
slower than $\sqrt{n}$, then the conditional loss or weighting estimators are
not $\sqrt{n}$ convergent and Wald-type confidence intervals are not
asymptotically valid [10]. Now we propose a doubly robust estimator that is
more robust to model misspecification (i.e., is consistent under milder
conditions) and allows for valid asymptotic inference even when data-adaptive
models that are not $\sqrt{n}$ convergent are used for nuisance parameter
estimation.

The doubly robust estimator combines modeling the probability of being in the
source data and the model for the conditional loss. For estimating $\psi$, the
doubly robust estimator is given by

$\displaystyle\widehat{\psi}_{DR}$
$\displaystyle=\frac{1}{n_{0}}\sum_{i=1}^{n}\Big{[}I(D_{i}=0)\widehat{h}(X_{i})+\frac{(1-\widehat{p}(X_{i}))I(D_{i}=1)}{\widehat{p}(X_{i})}(L(Y_{i},g(X^{*}_{i}))-\widehat{h}(X_{i}))\Big{]}.$
(7)

In Appendix A.1 we show how the doubly robust estimator is derived from the
efficient influence function under a non-parametric model for the data [16,
17]. Of note, the estimators presented in the previous section are special
cases of the doubly robust estimator: we obtain the weighting estimator by
setting $\widehat{h}(X_{i})=0$ identically for all $i$; we obtain the
conditional loss estimator by setting $\widehat{p}(X_{i})=1$ for all $i$.

###  3.3 Asymptotic properties of the doubly robust estimator

The presence of the estimators of nuisance parameters, $\widehat{p}(X)$ and
$\widehat{h}(X)$, in the definition of $\widehat{\psi}_{DR}$ complicates
derivation of large sample properties as standard central limit theorems and
laws of large numbers are not directly applicable. To derive the large sample
properties of the doubly robust estimator we define

$\displaystyle
H(p^{\prime}(X),h^{\prime}(X),\tau)=\tau^{-1}\left[I(D=0)h^{\prime}(X)+\frac{(1-p^{\prime}(X))I(D=1)}{p^{\prime}(X)}\big{(}L(Y,g(X^{*}))-h^{\prime}(X)\big{)}\right],$

for arbitrary functions $p^{\prime}(X),h^{\prime}(X),$ and $\tau$. The doubly
robust estimator can be written as
$\widehat{\psi}_{DR}=\frac{1}{n}\sum_{i=1}^{n}H(\widehat{p}(X_{i}),\widehat{h}(X_{i}),\widehat{\tau})$,
where $\widehat{p}(X_{i})$, and $\widehat{h}(X_{i})$ are as defined before,
and $\widehat{\tau}=\frac{1}{n}\sum_{j=1}^{n}I(D_{j}=0)$. Denote the limits of
$\widehat{p}(X)$ and $\widehat{h}(X)$ as $p^{*}(X)$ and $h^{*}(X)$,
respectively, and let $\tau_{0}=\Pr[D=0]$. Under correct model specification,
the limits are equal to $p^{*}(X)=\Pr[D=1|X]$ and
$h^{*}(X)=\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$.

We make the following four assumptions:

  * B1.

$H(\widehat{p}(X),\widehat{h}(X),\widehat{\tau})$ and its limit
$H(p^{*}(X),h^{*}(X),\tau_{0})$ fall in a Donsker class [18].

  * B2.

$||H(\widehat{p}(X),\widehat{h}(X),\widehat{\tau})-H(p^{*}(X),h^{*}(X),\tau_{0})||\overset{P}{\longrightarrow}0$.

  * B3.

(Finite second moment)
$\operatorname{\textnormal{\mbox{E}}}[H(p^{*}(X),h^{*}(X),\tau_{0})^{2}]<\infty$.

  * B4.

(Model double robustness) At least one of the models $\widehat{p}(X)$ or
$\widehat{h}(X)$ is correctly specified. That is, at least one of
$p^{*}(X)=\Pr[D=1|X]$ or
$h^{*}(X)=\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$ holds,
but not necessarily both.

Assumption B4 reflects the model double robustness property (i.e., it is
enough to correctly specify one of the models $\widehat{p}(X)$ and
$\widehat{h}(X)$, not both). If $\widehat{p}(X)$, $\widehat{h}(X)$,
$p^{*}(X)$, and $h^{*}(X)$ are Donsker, $\widehat{h}(X)$ and $h^{*}(X)$ are
uniformly bounded, and $\widehat{p}(X)$ and $p^{*}(X)$ are uniformly bounded
away from zero, then Assumption B1 holds by Donsker preservation theorems
[19]. The Donsker requirement restricts the complexity of the models
$\widehat{p}(X)$ and $\widehat{h}(X)$ and holds for many commonly used models
such as generalized linear models. To relax that assumption, sample splitting
can be used to allow for flexible machine learning algorithms to be used when
fitting $\widehat{h}(X)$ and $\widehat{p}(X)$ without violating Assumption B1
[10, 20, 21]. Using Assumptions B1 through B4, we prove the following theorem
in the Supplementary Web Appendix:

######  Theorem 3.

If assumptions B1 through B4 hold, then

  * 1.

(Consistency) $\widehat{\psi}_{DR}\overset{P}{\longrightarrow}\psi$.

  * 2.

$\widehat{\psi}_{DR}$ has the asymptotic representation

$\displaystyle\sqrt{n}(\widehat{\psi}_{DR}-\psi)$
$\displaystyle=\sqrt{n}\left(\frac{1}{n}\sum_{i=1}^{n}H(p^{*}(X_{i}),h^{*}(X_{i}),\tau_{0})-\operatorname{\textnormal{\mbox{E}}}\left[H(p^{*}(X),h^{*}(X),\tau_{0})\right]\right)+Re+o_{P}(1),$

where

$Re\leq\sqrt{n}\
O_{P}\bigg{(}||\widehat{h}(X)-\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]||_{2}^{2}\times||\widehat{p}(X)-\Pr[D=1|X]||_{2}^{2}\bigg{)}.$
(8)

Theorem 3 gives useful insights into the behaviour of the doubly robust
estimator and how the models used for estimation of $\Pr[D=1|X]$ and
$\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$ impact the large
sample properties of $\widehat{\psi}_{DR}$. Part 1 of Theorem 3 shows that
$\widehat{\psi}_{DR}$ is consistent if at least one of the models
$\widehat{p}(X)$ and $\widehat{h}(X)$ is correctly specified (model double
robustness [22]). In contrast, consistency of $\widehat{\psi}_{CL}$ relies on
the estimator $\widehat{h}(X)$ to be correctly specified and consistency of
$\widehat{\psi}_{IW}$ relies on the estimator $\widehat{p}(X)$ to be correctly
specified. Furthermore, the first term on the right hand side of the
asymptotic representation of $\widehat{\psi}_{DR}$ (in Part 2 of Theorem 3)
has the estimators $\widehat{p}(X)$ and $\widehat{h}(X)$ replaced by their
limits; thus, by the Central Limit Theorem, that term converges asymptotically
to a normal distribution at $\sqrt{n}$ rate.

The inequality displayed in (8) shows that $\widehat{\psi}_{DR}$ converges at
$\sqrt{n}$-rate if the combined rate of convergence of $\widehat{p}(X)$ and
$\widehat{h}(X)$ is at least $\sqrt{n}$ (rate double robustness [22]).
Contrast this to $\widehat{\psi}_{CL}$ and $\widehat{\psi}_{IW}$ that inherit
the rate of convergence of $\widehat{h}(X)$ and $\widehat{p}(X)$, respectively
[23]. This implies that estimation procedures such as generalized additive
models and the highly adaptive Lasso [24] can be used to estimate
$\widehat{h}(X)$ and $\widehat{p}(X)$ in $\widehat{\psi}_{DR}$ while still
being $\sqrt{n}$ convergent (both have rate of convergence slower than
$\sqrt{n}$ but faster than $n^{1/4}$ [25, 24]). For inference, the variance of
$\widehat{\psi}_{DR}$ can be estimated using either a sandwich variance
estimator [26] or using the bootstrap.

##  4 Simulation study

###  4.1 Setup

We used simulations to compare the performance of the estimators for the
target population risk $\psi$: the weighting estimator $\widehat{\psi}_{IW}$,
the conditional loss estimator $\widehat{\psi}_{CL}$, the doubly robust
estimator $\widehat{\psi}_{DR}$, and a naïve empirical estimator that
estimates the target population risk using the loss calculated using only the
source data
$\widehat{\psi}_{S}=\frac{1}{n_{1}}\sum_{i=1}^{n}I(D_{i}=1)L(Y_{i},g(X^{*}_{i}))$.

We simulated the covariate vector $X$ from a $10$ dimensional mean zero
multivariate normal distribution with element $(i,j)$ of the covariance matrix
equal to $0.5^{|i-j|}$. We simulated the probability of being from the source
data from a Bernoulli distribution with parameter
$\Pr[D=1|X]=\frac{exp(-0.3+0.2\sum_{i=1}^{3}X^{(i)}+0.3\sum_{i=1}^{3}(X^{(i)})^{2})}{1+exp(-0.3+0.2\sum_{i=1}^{3}X^{(i)}+0.3\sum_{i=1}^{3}(X^{(i)})^{2})}$,
where $X^{(i)}$ is the $i$th component of the covariate vector $X$. This
resulted in $\Pr[D=1]=0.61$. We simulated the outcome from a Bernoulli
distribution with parameter
$\Pr[Y=1|X]=\frac{exp(-0.3+0.2\sum_{i=1}^{3}X^{(i)}+0.3\sum_{i=1}^{3}(X^{(i)})^{2})}{1+exp(-0.3+0.2\sum_{i=1}^{3}X^{(i)}+0.3\sum_{i=1}^{3}(X^{(i)})^{2})}$.
The source data was split into a training set (2/3 of the source data) and a
test set consisting of the remaining observations from the sample from the
source population. The prediction model $g(X^{*})$ that we assessed the
performance of was a main effects logistic regression model fit using the
training data from the sample from the source population. No data from the
target population was used to estimate $g(X^{*})$, so all observations in the
target data and the test set from the sample from the source population were
used to estimate the Brier loss in the target population (the evaluation
measure we focus on). For each simulation we sampled a total of $1000$
observations (both $D=1$ and $D=0$), and to compare the four estimators we
conducted $1000$ simulations.

For binary outcomes, expanding the square shows that to estimate
$h(X)=\operatorname{\textnormal{\mbox{E}}}[(Y-g(X^{*}))^{2}|X,D=1]$ it is
enough to estimate $\Pr[Y=1|X,D=1]$, and this is the approach we used in the
simulations. To evaluate the impact of model specification on performance, we
compared the four estimators using different combinations of correctly
specified and misspecified models for $\Pr[Y=1|X,D=1]$ and $\Pr[D=1|X]$. For
the correctly specified models, we estimated $\Pr[Y=1|X,D=1]$ and $\Pr[D=1|X]$
using a main effects logistic regression model that included both linear and
quadratic main effects of all covariates. To misspecify the models for
$\Pr[Y|X,D=1]$ and $\Pr[D=1|X]$ we fit a logistic regression model that only
included linear main effects (i.e., both models failed to include the
quadratic main effects). We also implemented a doubly robust estimator that
used generalized additive models to estimate both $\Pr[Y=1|X,D=1]$ and
$\Pr[D=1|X]$. We fit the generalized additive models using the mgcv package in
R [27] entering all covariates as splines using the default options in the gam
function.

Based on the theoretical properties of the estimators, we expected the naïve
empirical estimator to be biased. We also expected the weighting estimator to
be biased when the model for source study participation was misspecified, and
the conditional loss estimator to be biased when the model for the conditional
loss was incorrectly specified. Finally, we expected the doubly robust
estimator to be unbiased when at least one of the models was correctly
specified.

###  4.2 Simulation results

Figure 1 shows boxplots of the results (exact numerical values are presented
in Table S1 in Web Appendix). As expected, the naïve empirical estimator was
biased downwards with a relative bias of $12\%$. Under misspecification of
$\widehat{p}(X)$, the weighting estimator was biased but the doubly robust
estimator is unbiased. Under misspecification of $\widehat{h}(X)$, the
conditional loss estimator was biased, but the doubly robust estimator was
unbiased. When both models $\widehat{p}(X)$ and $\widehat{h}(X)$ were
misspecified all estimators, including the doubly robust estimator, were
biased. When a generalized additive model was used to estimate both
$\Pr[Y=1|X,D=1]$ and $\Pr[D=1|X]$, the doubly robust estimator was
approximately unbiased (relative bias of 1.1%).

##  5 Assessing the performance of a model for lung cancer diagnosis

We applied the proposed methods to evaluate the performance of a prediction
model for lung cancer diagnosis in a target population of people eligible for
lung cancer screening in the US.

###  5.1 Study design and data

We used source population data from the National Lung Screening Trial (NLST),
a large clinical trial comparing the effect of screening for lung cancer using
computerized tomography (CT) versus chest X-ray [28]. The NLST enrolled
participants from 2002-2004 and followed them until the end of December 2009.
The trial found that CT screening reduced lung cancer specific mortality and
it has had a large influence on guidelines for lung cancer screening [29, 30].
We used data from the CT arm of the NLST ($n_{1}=25,825$ after removing 897
observations with missing data). We used lung cancer diagnosis within six
years as the outcome.

We defined the target population as the population of non-institutionalized
U.S. adults who met the eligibility criteria of the NLST. To obtain
information on the covariate distribution of the target population, we used
data from the 2009-2010 cycle of NHANES, which is sampled to be representative
of the non-institutionalized US population. We focused on the subset of NHANES
participants that met the NLST eligibility criteria (55–74 years old with no
history of lung cancer who currently smoke or have quit within the past 15
years and that have greater than 30 pack-year history) and were included in a
smoking sub-study (and therefore had information on smoking habits needed to
assess eligibility). This resulted in a sample size of $n_{0}=189$ after
removing 4 observations with missing data. The NLST eligibility criteria are
almost identical to the criteria used for recommending people for lung cancer
screening in the 2014 United States Preventive Services Task Force (USPSTF)
guidelines [29] and very similar to the updated criteria in the 2021 USPSTF
guidelines for lung cancer screening [30]. Because NHANES is a cross-sectional
study, it does not collect long-term follow-up information on lung cancer
diagnosis; that is to say, no outcome data are available in the sample from
the target population. Therefore, we cannot build a prediction model and/or
evaluate the performance of a prediction model using only data from NHANES.

Table S2 in Web Appendix C summarizes variables from the NSLT and NHANES
datasets that we used in our analyses. The data suggest the presence of
substantial differences between the population underlying the NLST trial and
the target population: participants in the NLST were less ethnically and
racially diverse and are more educated compared to the national population of
people eligible for NLST (represented by the subset of the NHANES data that
met the NLST eligibility criteria). NHANES used a multi-cluster sampling
design where each observation is associated with a sampling weight (accounting
for oversampling of certain subgroups, survey non-response, and post-
stratification adjustments). In subsection 5.2 we show how the conditional
loss, weighting, and doubly robust estimators can be modified to incorporate
the survey sampling weights and account for the multi-stage clustering.

###  5.2 Accounting for the complex survey design of the NHANES

#### Incorporating the NHANES survey sampling weights:

Let $w_{i}$ be the sampling weight for observation $i$ (where all observations
in the sample from the source population get a weight of $1$). To account for
the sampling weights, we modify the conditional loss estimator as follows:

$\widehat{\psi}_{CL}^{w}=\frac{1}{\sum_{i=1}^{n}w_{i}I(D_{i}=0)}\sum_{i=1}^{n}w_{i}I(D_{i}=0)\widehat{h}(X_{i}).$
(9)

The estimator $\widehat{h}(X)$ for
$\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$ is obtained using
only the sample from the source population; therefore, it does not need to
incorporate the survey sampling weights.

Furthermore, we modify the weighting estimator as follows:

$\widehat{\psi}_{IW}^{w}=\frac{1}{\sum_{i=1}^{n}w_{i}I(D_{i}=0)}\sum_{i=1}^{n}\frac{I(D_{i}=1)(1-\widehat{\tilde{p}}(X_{i}))}{\widehat{\tilde{p}}(X_{i})}L(Y_{i},g(X^{*}_{i})),$
(10)

where $\widehat{\tilde{p}}(X)$ is an estimator for $\Pr[D=1|X]$ that accounts
for the sampling weights in the target population (e.g., using a weighted
logistic regression model). Last, we modify the doubly robust estimator as
follows:

$\displaystyle\widehat{\psi}_{DR}^{w}$
$\displaystyle=\frac{\sum_{i=1}^{n}\Big{[}w_{i}I(D_{i}=0)\widehat{h}(X_{i})+\frac{(1-\widehat{\tilde{p}}(X_{i}))I(D_{i}=1)}{\widehat{\tilde{p}}(X_{i})}(L(Y_{i},g(X^{*}_{i}))-\widehat{h}(X_{i}))\Big{]}}{\sum_{i=1}^{n}w_{i}I(D_{i}=0)}.$
(11)

#### Model specification:

We evaluated the performance a logistic regression model with the following
predictors entered as main effects: age, BMI, race (Black, White, Hispanic,
other), education (less than high school, high school graduate, associate’s
degree/some college, college graduate), personal history of cancer, smoking
status, smoking intensity, duration of smoking, and smoking quit time. We
split the NLST data into a training and a test set of equal size and the model
was fit on the training set and the Brier score [5] in the target population
was estimated using the test set from the NLST data and all the NHANES data.
To implement the conditional loss, weighting, and doubly robust estimators we
used a main effects logistic regression models to estimate $\Pr[D=1|X]$
weighted by the survey sampling weights; to estimate
$\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$ we used a main
effects logistic regression model.

#### Inference:

To quantify uncertainty about risk estimates in the target population, we used
the non-parametric bootstrap with $1000$ bootstrap samples. To account for the
multi-stage clustering design, we resampled the data consistent with the
NHANES sampling design [31]. To obtain the standard deviation of the naïve
empirical estimator, we used the non-parametric bootstrap with $1000$
bootstrap samples.

###  5.3 Results

Table 1 shows estimates and associated standard errors for the Brier risk in
the target population from the naïve empirical, conditional loss, weighting,
and doubly robust estimators. The conditional loss, weighting, and doubly
robust estimators produced similar estimates between them, but substantially
greater than those of the naïve empirical estimator. The weighting estimator
had greater standard errors than the doubly robust estimator (by 12%) as well
as the conditional loss estimator (by 53%).

In Web Appendix D we present the results from an additional proof-of-principle
analysis where the NLST data were artificially split into a sample from the
source population and a sample from the target population (under the covariate
shift assumption). In that analysis, by construction, outcome information was
available from the target population sample, allowing us to compare the
proposed estimators to an “oracle” estimator that used outcome information
from the target population. The results show similar trends as in the
simulations presented in Section 4, the naïve empirical estimator produced
estimates that differed from those of the oracle estimator while the three
estimators of the target population Brier risk that account for covariate
shift produced estimates similar to the oracle estimator. Furthermore, the
weighting estimator had greater variance than the doubly robust and
conditional loss estimators.

##  6 Discussion

If prediction error modifiers have a different distribution between the source
population from which data for model development are obtained, and the target
population where the model will be applied, then measures of model performance
calculated using data only from the source population may not directly apply
to the target population [6]. Hence, methods that can adjust for between-
population differences in prediction error modifiers are needed when assessing
model performance in the target population.

We proposed and developed theoretical properties of a doubly robust estimator
for the risk in the target population that is more robust to model
misspecification than the previously proposed weighting estimator and can be
used with data-adaptive estimators of nuisance parameters that converge at a
rate slower than $\sqrt{n}$. As a special case of the doubly robust estimator,
we also proposed a novel conditional loss estimator. We compared the doubly
robust and conditional loss estimators to the weighting estimator using
simulations and data on lung cancer screening and found that they performed
favorable compared to previously proposed inverse weighting estimators.

In our setup, the source data is split into a training set used to build the
prediction model and a test set that is combined with the sample from the
target population to estimate the performance of the prediction model.
Alternatively, we may be interested in assessing model performance of an
already established prediction model (built using an external data source). In
that case, the whole sample from the source population (not just the test set)
can be used for assessment of model performance and the methods developed in
this manuscript are still valid.

## References

  * [1] DR Cox, Christiana Kartsonaki, and Ruth H Keogh.  Big data: Some statistical issues.  Statistics & Probability Letters, 136:111–115, 2018.
  * [2] Hidetoshi Shimodaira.  Improving predictive inference under covariate shift by weighting the log-likelihood function.  Journal of Statistical Planning and Inference, 90(2):227–244, 2000\.
  * [3] Steffen Bickel, Michael Brückner, and Tobias Scheffer.  Discriminative learning under covariate shift.  Journal of Machine Learning Research, 10(Sep):2137–2155, 2009.
  * [4] Jiayuan Huang, Arthur Gretton, Karsten Borgwardt, Bernhard Schölkopf, and Alex J Smola.  Correcting sample selection bias by unlabeled data.  In Advances in Neural Information Processing Systems, pages 601–608, 2007.
  * [5] Glenn W Brier.  Verification of forecasts expressed in terms of probability.  Monthly Weather Review, 78(1):1–3, 1950.
  * [6] Jon A. Steingrimsson, Constantine Gatsonis, and Issa J. Dahabreh.  Transporting a prediction model for use in a new target population.  arXiv preprint arXiv:2101.11182, 2021.
  * [7] Bianca Zadrozny.  Learning and evaluating classifiers under sample selection bias.  In Proceedings of the Twenty-First International Conference on Machine Learning, page 114. ACM, 2004.
  * [8] Masashi Sugiyama, Matthias Krauledat, and Klaus-Robert MÃžller.  Covariate shift adaptation by importance weighted cross validation.  Journal of Machine Learning Research, 8(May):985–1005, 2007.
  * [9] James M Robins, Andrea Rotnitzky, and Lue Ping Zhao.  Estimation of regression coefficients when some regressors are not always observed.  Journal of the American Statistical Association, 89(427):846–866, 1994.
  * [10] Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins.  Double/debiased machine learning for treatment and structural parameters.  The Econometrics Journal, 21(1):C1–C68, 2018.
  * [11] Issa J Dahabreh, Sebastien JP A Haneuse, James M Robins, Sarah E Robertson, Ashley L Buchanan, Elizabeth A Stuart, and Miguel A Hernán.  Study designs for extending causal inferences from a randomized trial to a target population.  American journal of epidemiology, 190(8):1632–1642, 2021.
  * [12] Jon A Steingrimsson.  Extending prediction models for use in a new target population with failure time outcomes.  Biostatistics, 2022.
  * [13] James M Robins and Andrea Rotnitzky.  Recovery of information and adjustment for dependent censoring using surrogate markers.  In AIDS Epidemiology, pages 297–331. Springer, 1992.
  * [14] Daniel G Horvitz and Donovan J Thompson.  A generalization of sampling without replacement from a finite universe.  Journal of the American Statistical Association, 47(260):663–685, 1952.
  * [15] Jerome Friedman, Trevor Hastie, Robert Tibshirani, et al.  The elements of statistical learning, volume 1.  Springer series in statistics New York, 2001.
  * [16] Aaron Fisher and Edward H Kennedy.  Visually communicating and teaching intuition for influence functions.  The American Statistician, pages 1–11, 2020.
  * [17] Mark J Van der Laan, MJ Laan, and James M Robins.  Unified methods for censored longitudinal data and causality.  Springer Science & Business Media, 2003.
  * [18] Richard M Dudley.  Uniform central limit theorems, volume 142.  Cambridge University Press, 2014.
  * [19] Michael R Kosorok.  Introduction to empirical processes and semiparametric inference.  Springer Science & Business Media, 2007.
  * [20] Peter J Bickel and Yaacov Ritov.  Estimating integrated squared density derivatives: sharp best order of convergence estimates.  Sankhyā: The Indian Journal of Statistics, Series A, pages 381–393, 1988.
  * [21] James Robins, Lingling Li, Eric Tchetgen Tchetgen, and Aad van der Vaart.  Higher order influence functions and minimax estimation of nonlinear functionals.  In Probability and statistics: essays in honor of David A. Freedman, pages 335–421. Institute of Mathematical Statistics, 2008.
  * [22] Andrea Rotnitzky, Ezequiel Smucler, and James M Robins.  Characterization of parameters with a mixed bias property.  Biometrika, 108(1):231–238, 2021.
  * [23] Iván Díaz.  Machine learning in the estimation of causal effects: targeted minimum loss-based estimation and double/debiased machine learning.  Biostatistics, 21(2):353–358, 2020.
  * [24] David Benkeser and Mark Van Der Laan.  The highly adaptive lasso estimator.  In 2016 IEEE International Conference on Data Science and Advanced Analytics (DSAA), pages 689–696. IEEE, 2016.
  * [25] Joel L Horowitz.  Semiparametric and nonparametric methods in econometrics, volume 12.  Springer, 2009.
  * [26] Leonard A Stefanski and Dennis D Boos.  The calculus of m-estimation.  The American Statistician, 56(1):29–38, 2002.
  * [27] Simon Wood and Maintainer Simon Wood.  Package ‘mgcv’.  R package version, 1(29):729, 2015.
  * [28] National Lung Screening Trial Research Team.  Reduced lung-cancer mortality with low-dose computed tomographic screening.  New England Journal of Medicine, 365(5):395–409, 2011.
  * [29] Virginia A Moyer.  Screening for lung cancer: Us preventive services task force recommendation statement.  Annals of internal medicine, 160(5):330–338, 2014.
  * [30] Alex H Krist, Karina W Davidson, Carol M Mangione, Michael J Barry, Michael Cabana, Aaron B Caughey, Esa M Davis, Katrina E Donahue, Chyke A Doubeni, Martha Kubik, et al.  Screening for lung cancer: Us preventive services task force recommendation statement.  JAMA, 325(10):962–970, 2021.
  * [31] Jun Shao.  Impact of the bootstrap on sample surveys.  Statistical Science, 18(2):191–198, 2003.

## Tables and Figures

Table 1: Estimates for the Brier risk in the target population and the associated estimates of standard error estimated using data from the National Lung Screening Trial and the NHANES survey. | $\widehat{\psi}_{S}$ | $\widehat{\psi}_{IW}^{w}$ | $\widehat{\psi}_{CL}^{w}$ | $\widehat{\psi}_{DR}^{w}$
---|---|---|---|---
Estimator for Brier risk | 0.039 | 0.052 | 0.049 | 0.049
Standard error estimator | 0.0016 | 0.0075 | 0.0049 | 0.0067
$\widehat{\psi}_{S}$ is the naïve empirical estimator,
$\widehat{\psi}_{IW}^{w}$ is the weighting estimator,
$\widehat{\psi}_{CL}^{w}$ is the conditional loss estimator, and
$\widehat{\psi}_{DR}^{w}$ is the doubly robust estimator.

Figure 1: Boxplots of the difference between the estimated and the true Brier
loss (estimated numerically) for the naïve empirical, weighting (IW),
conditional loss (CL), and doubly robust (DR) estimators. Values closer to
zero indicate better performance and the simulation setup is described in
Section 4.  Here, $h(X)$ denotes a model for $\Pr[Y=1|X,D=1]$ and $p(X)$
denotes a model for $\Pr[D=1|X]$. Corr indicates correct specification of the
logistic regression models needed for implementation of the estimators. Miss
indicates a misspecified logistic regression model. GAM indicates that a
generalized additive model was used to estimate both $\Pr[Y=1|X,D=1]$ and
$\Pr[D=1|X]$.

##  Appendix A Proofs of key results

###  A.1 Efficient influence function and doubly robust estimator

Now we will show that the influence function under a non-parametric model for
the observable data is

$\displaystyle\chi_{p_{0}}^{1}$
$\displaystyle=\frac{I(D=0)}{\Pr[D=0]}\big{\\{}E[L(Y,g(X^{*}))|X,D=1]-\psi\big{\\}}$
$\displaystyle+\frac{\Pr[D=0|X]I(D=1)}{\Pr[D=0]\Pr[D=1|X]}\big{\\{}L(Y,g(X^{*}))-\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]\big{\\}}.$

As the influence function under the non-parametric model is unique it is also
the efficient influence function.

To prove that $\chi_{p_{0}}^{1}$ is the influence function we show that
$\chi_{p_{0}}^{1}$ is a solution to

$\frac{d}{dt}\psi_{p_{t}}\Bigr{|}_{t=0}=\operatorname{\textnormal{\mbox{E}}}_{p_{0}}[\chi_{p_{0}}^{1}g_{p_{0}}],$

where the left hand side of the above equation is the pathwise derivative of
the target parameter and $g_{p_{0}}$ is the score of the observable data [17].
Here, the subscript $p_{0}$ denotes that the expectation is calculated under
the true data law. For notational convenience, define
$h(X)=\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$ and
$O=(Y,X,D)$.

The pathwise derivative of $\psi_{p_{t}}$ with respect to t where the
subscript $p_{t}$ denotes the dependence of $\psi_{p_{t}}$ on a one-
dimensional parametric sub-model $p_{t}$ indexed by $t\in[0,1]$ is given by

$\displaystyle\frac{\partial}{\partial t}\psi_{p_{t}}\Bigr{|}_{t=0}$
$\displaystyle=\frac{d}{dt}\operatorname{\textnormal{\mbox{E}}}_{pt}[\operatorname{\textnormal{\mbox{E}}}_{pt}[L(Y,g(X^{*}))|X,D=1]|D=0]\Bigr{|}_{t=0}$
$\displaystyle=\frac{\partial}{\partial
t}\operatorname{\textnormal{\mbox{E}}}_{pt}[\operatorname{\textnormal{\mbox{E}}}_{p_{0}}[L(Y,g(X^{*}))|X,D=1]|D=0]\Bigr{|}_{t=0}$
$\displaystyle+\operatorname{\textnormal{\mbox{E}}}_{p_{0}}\left[\frac{\partial}{\partial
t}\operatorname{\textnormal{\mbox{E}}}_{pt}[L(Y,g(X^{*}))|X,D=1]\Bigr{|}_{t=0}\Bigr{|}D=0\right]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{p_{0}}\left[\big{\\{}h(X)-\psi\big{\\}}\frac{I(D=0)}{\Pr[D=0]}g_{Y,X,D}(O)\right]$
$\displaystyle+\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}_{p_{0}}\left[\big{\\{}L(Y,g(X^{*}))-h(X)\big{\\}}\frac{I(D=1)\Pr[D=0|X]}{\Pr[D=1|X]}g_{Y,X,D}(O)\right]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{p_{0}}\Bigg{[}\bigg{(}\big{\\{}h(X)-\psi\big{\\}}\frac{I(D=0)}{\Pr[D=0]}$
$\displaystyle+\frac{1}{\Pr[D=0]}\big{\\{}L(Y,g(X^{*}))-h(X)\big{\\}}\Big{\\{}\frac{I(D=1)\Pr[D=0|X]}{\Pr[D=1|X]}\Big{\\}}\bigg{)}g_{Y,X,D}(O)\Bigg{]}.$

It follows that the influence function is

$\displaystyle\chi_{p_{0}}^{1}$
$\displaystyle=\frac{I(D=0)}{\Pr[D=0]}\big{\\{}E_{p_{0}}[L(Y,g(X^{*}))|X,D=1]-\psi\big{\\}}$
$\displaystyle+\frac{\Pr[D=0|X]I(D=1)}{\Pr[D=0]\Pr[D=1|X]}\big{\\{}L(Y,g(X^{*}))-\operatorname{\textnormal{\mbox{E}}}_{p_{0}}[L(Y,g(X^{*}))|X,D=1]\big{\\}}$

###  A.2 Proofs of asymptotic properties of the doubly robust estimator

#### Consistency of $\widehat{\psi}_{DR}.$

Define the limits of $\widehat{p}(X)$ and $\widehat{h}(X)$ (assumed to exist)
as $p^{*}(X)$ and $h^{*}(X)$, respectively. Under correct model specification
the limits are equal to $p^{*}(X)=\Pr[D=1|X]$ and
$h^{*}(X)=\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$.

The doubly robust estimator $\widehat{\psi}_{DR}$ converges in probability to

$\displaystyle\widehat{\psi}_{DR}\overset{P}{\longrightarrow}\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}\left[I(D=0)h^{*}(X)+\frac{(1-p^{*}(X))I(D=1)}{p^{*}(X)}(L(Y,g(X^{*}))-h^{*}(X))\right]$

Now we show that the quantity on the right-hand side is equal to $\psi$ under
assumptions B1- B4. First consider the case where $\widehat{p}(X)$ is
correctly specified, that is $p^{*}(X)=\Pr[D=1|X]$ (we do not make the
assumption that the limit $h^{*}(X)$ is equal to
$\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$).

Algebra shows that
$\psi=\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}\big{[}\frac{\Pr[D=0|X]I(D=1)}{\Pr[D=1|X]}L(Y,g(X^{*}))\big{]}$.
Using that we get

$\displaystyle\widehat{\psi}_{DR}$
$\displaystyle\overset{P}{\longrightarrow}\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}\left[I(D=0)h^{*}(X)+\frac{\Pr[D=0|X]I(D=1)}{\Pr[D=1|X]}(L(Y,g(X^{*}))-h^{*}(X))\right]$
$\displaystyle=\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}\left[I(D=0)h^{*}(X)-\frac{\Pr[D=0|X]I(D=1)}{\Pr[D=1|X]}h^{*}(X)\right]+\psi$
$\displaystyle=\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}\left[\operatorname{\textnormal{\mbox{E}}}\left[I(D=0)h^{*}(X)-\frac{\Pr[D=0|X]I(D=1)}{\Pr[D=1|X]}h^{*}(X)\bigg{|}X\right]\right]+\psi$
$\displaystyle=\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}\left[h^{*}(X)\operatorname{\textnormal{\mbox{E}}}[I(D=0)|X]-\frac{\Pr[D=0|X]}{\Pr[D=1|X]}h^{*}(X)\operatorname{\textnormal{\mbox{E}}}[I(D=1)|X]\right]+\psi$
$\displaystyle=\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}\left[h^{*}(X)\Pr[D=0|X]-\frac{\Pr[D=0|X]}{\Pr[D=1|X]}h^{*}(X)\Pr[D=1|X]\right]+\psi$
$\displaystyle=\psi.$

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

$h^{*}(X)=\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$

and we do not make the assumptions that the limit $p^{*}(X)$ is equal to
$\Pr[D=1|X]$. Algebra shows

$\psi=\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}\big{[}I(D=0)\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]\big{]}.$

Hence,

$\displaystyle\widehat{\psi}_{DR}$
$\displaystyle\overset{P}{\longrightarrow}\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}\Bigg{[}I(D=0)\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$
$\displaystyle\quad\quad\quad+\frac{(1-p^{*}(X))I(D=1)}{p^{*}(X)}(L(Y,g(X^{*}))-\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1])\Bigg{]}$
$\displaystyle=\psi$

#### Asymptotic representation of $\widehat{\psi}_{DR}$.

For a random variable $W$ define

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

Using this notation, we rewrite the asymptotic representation as

$\displaystyle\sqrt{n}(\widehat{\psi}_{DR}-\psi)$
$\displaystyle=\mathbb{G}_{n}(H(\widehat{p}(X),\widehat{h}(X),\widehat{\tau}))-\mathbb{G}_{n}(H(p^{*}(X),h^{*}(X),\tau_{0}))$
(S1)
$\displaystyle\quad\quad\quad+\mathbb{G}_{n}(H(p^{*}(X),h^{*}(X),\tau_{0}))$
(S2)
$\displaystyle\quad\quad\quad+\sqrt{n}\left(\operatorname{\textnormal{\mbox{E}}}[H(\widehat{p}(X),\widehat{h}(X),\widehat{\tau})]-\psi\right)$
(S3)

For term S1 in the above equation, the Donsker assumption (Assumption A1)
implies that [18]

$\mathbb{G}_{n}(H(\widehat{p}(X),\widehat{h}(X),\widehat{\tau}))-\mathbb{G}_{n}(H(p^{*}(X),h^{*}(X),\tau_{0}))=o_{P}(1).$

By defining

$Re=\sqrt{n}(\operatorname{\textnormal{\mbox{E}}}[H(\widehat{p}(X),\widehat{h}(X),\widehat{\tau})]-\psi),$

we have

$\sqrt{n}(\widehat{\psi}_{DR}-\psi)=\sqrt{n}\left(\frac{1}{n}\sum_{i=1}^{n}\big{(}H(p^{*}(X_{i}),h^{*}(X_{i}),\tau_{0})-\operatorname{\textnormal{\mbox{E}}}\left[H(p^{*}(X),h^{*}(X),\tau_{0})\right]\big{)}\right)+Re+o_{P}(1).$

Now we calculate the upper bound of $Re$. First rewrite

$\displaystyle
n^{-1/2}Re=\frac{1}{\Pr[D=0]}\bigg{(}\underbrace{\operatorname{\textnormal{\mbox{E}}}\left[I(D=0)\widehat{h}(X)\right]}_{R_{1}}+\underbrace{\operatorname{\textnormal{\mbox{E}}}\left[\frac{(1-\widehat{p}(X))I(D=1)}{\widehat{p}(X)}\\{L(Y,g(X^{*}))-\widehat{h}(X)\\}\right]}_{R_{2}}\bigg{)}-\psi.$

The first term $R_{1}$ can be rewritten as:

$\displaystyle R_{1}$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}\left[I(D=0)\widehat{h}(X)\right]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}\left[\operatorname{\textnormal{\mbox{E}}}[I(D=0)\widehat{h}(X)|X]\right]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}\left[\operatorname{\textnormal{\mbox{E}}}[I(D=0)|X]\widehat{h}(X)\right]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}[\Pr[D=0|X]\widehat{h}(X)].$

We rewrite term $R_{2}$ as:

$\displaystyle R_{2}$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}\left[\frac{(1-\widehat{p}(X))I(D=1)}{\widehat{p}(X)}\\{L(Y,g(X^{*}))-\widehat{h}(X)\\}\right]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}\left[\operatorname{\textnormal{\mbox{E}}}\left[\frac{(1-\widehat{p}(X))I(D=1)}{\widehat{p}(X)}\\{L(Y,g(X^{*}))-\widehat{h}(X)\\}\bigg{|}X\right]\right]\
\ \ \ $
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}\left[\frac{(1-\widehat{p}(X))}{\widehat{p}(X)}\operatorname{\textnormal{\mbox{E}}}\left[\frac{I(D=1)}{\Pr[D=1|X]}\Pr[D=1|X]\\{L(Y,g(X^{*}))-\widehat{h}(X)\\}\bigg{|}X\right]\right]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}\left[\frac{(1-\widehat{p}(X))}{\widehat{p}(X)}\operatorname{\textnormal{\mbox{E}}}\left[\Pr[D=1|X]\\{L(Y,g(X^{*}))-\widehat{h}(X)\\}|X,D=1\right]\right]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}\left[\frac{(1-\widehat{p}(X))}{\widehat{p}(X)}\Pr[D=1|X]\\{\operatorname{\textnormal{\mbox{E}}}\left[L(Y,g(X^{*}))|X,D=1\right]-\widehat{h}(X)\\}\right]$

Finally, we rewrite,

$\displaystyle\psi$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}[\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]|D=0]$
$\displaystyle=\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}[I(D=0)\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]]$
$\displaystyle=\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}[\operatorname{\textnormal{\mbox{E}}}[I(D=0)\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]|X]]$
$\displaystyle=\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}[\Pr[D=0|X]\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]]$

Combining the above gives

$\displaystyle n^{-1/2}Re$
$\displaystyle=\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}[\Pr[D=0|X]\widehat{h}(X)]$
$\displaystyle+\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}\left[\frac{(1-\widehat{p}(X))}{\widehat{p}(X)}\Pr[D=1|X]\\{\operatorname{\textnormal{\mbox{E}}}\left[L(Y,g(X^{*}))|X,D=1\right]-\widehat{h}(X)\\}\right]$
$\displaystyle-\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}[\Pr[D=0|X]\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]]$
$\displaystyle=\frac{1}{\Pr[D=0]}\operatorname{\textnormal{\mbox{E}}}\big{[}-\Pr[D=0|X]\\{\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]-\widehat{h}(X)\\}$
$\displaystyle+\frac{(1-\widehat{p}(X))}{\widehat{p}(X)}\Pr[D=1|X]\\{\operatorname{\textnormal{\mbox{E}}}\big{[}L(Y,g(X^{*}))|X,D=1\big{]}-\widehat{h}(X)\\}\big{]}$
$\displaystyle=\frac{1}{\Pr[D=0]}E\Bigg{[}\\{\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]-\widehat{h}(X)\\}$
$\displaystyle\times\left\\{\frac{(1-\widehat{p}(X))}{\widehat{p}(X)}\Pr[D=1|X]-\Pr[D=0|X]\right\\}\Bigg{]}$

Using the Cauchy-Schwartz inequality we get.

$\displaystyle Re$
$\displaystyle\leq\sqrt{n}\frac{1}{\Pr[D=0]}\left(\operatorname{\textnormal{\mbox{E}}}[\\{\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]-\widehat{h}(X)\\}^{2}]\right)^{1/2}$
$\displaystyle\times\left(\operatorname{\textnormal{\mbox{E}}}\left[\left\\{\frac{(1-\widehat{p}(X))}{\widehat{p}(X)}\Pr[D=1|X]-\Pr[D=0|X]\right\\}^{2}\right]\right)^{1/2}$
$\displaystyle\leq$ $\displaystyle\sqrt{n}\
O_{P}\bigg{(}||\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]-\widehat{h}(X)||_{2}^{2}\times||\widehat{p}(X)-\Pr[D=1|X]||_{2}^{2}\bigg{)}.$

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

$\displaystyle\sqrt{n}(\widehat{\psi}_{DR}-\psi)$
$\displaystyle=\sqrt{n}\Bigg{(}\frac{1}{n}\sum_{i=1}^{n}H(\Pr[D=1|X_{i}],\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|D=1,X_{i}],\tau_{0})$
$\displaystyle-\operatorname{\textnormal{\mbox{E}}}[H(\Pr[D=1|X],\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|D=1,X],\tau_{0})]\Bigg{)}+o_{P}(1),$

By the central limit theorem,

$\sqrt{n}\left(\frac{1}{n}\sum_{i=1}^{n}H(p^{*}(X_{i}),h^{*}(X_{i}),\tau_{0})-\operatorname{\textnormal{\mbox{E}}}\left[H(p^{*}(X_{i}),h^{*}(X_{i}),\tau_{0})\right]\right)\overset{d}{\longrightarrow}N(0,Var[H(p^{*}(X),h^{*}(X),\tau_{0})])$

completing the proof.

##  Appendix B Additional simulation results

Table S1: Numerical results from the simulations described in Section 4 of the main text.  Estimator | Average of estimates | $\sqrt{n}\times Bias$ | $\sqrt{n}\times SD$ | Relative Bias
---|---|---|---|---
Naïve empirical | 0.2279 | -0.577 | 0.225 | -12%
W Corr | 0.2597 | 0.00316 | 0.285 | 0.067%
W Miss | 0.2520 | -0.136 | 0.264 | -2.9%
CL Corr | 0.2597 | 0.00379 | 0.256 | 0.080%
CL Miss | 0.2351 | -0.445 | 0.207 | -9.4%
DR Corr | 0.2595 | 0.000257 | 0.269 | 0.0054%
DR Miss $\widehat{p}(X)$ | 0.2596 | 0.00257 | 0.259 | 0.054%
DR Miss $\widehat{h}(X)$ | 0.2591 | -0.00733 | 0.268 | -0.15%
DR Miss Both | 0.2430 | -0.302 | 0.239 | -6.4%
DR GAM | 0.2565 | -0.05030 | 0.299 | -1.1%
Average of estimates, estimated bias, estimated standard deviation (SD), and
estimated relative bias for the naïve empirical, weighting (W), conditional
loss (CL), and doubly robust (DR) estimators. The simulation setup is
described in Section 4 of the main paper. $\sqrt{n}$ is the number of
observations used to fit the estimators for the target population risk (test
set of $D=1$ and $D=0$). Here, $\widehat{h}(X)$ is a model for
$\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$ and
$\widehat{p}(X)$ denotes a model for $\Pr[D=1|X]$. Relative bias is calculated
as $\frac{estimator-truth}{truth}$. Corr indicates correct specification of
all models needed for implementation of the estimator. Miss indicates a
misspecified model. DR GAM indicates that a generalized additive model was
used to estimate both
$\operatorname{\textnormal{\mbox{E}}}[L(Y,g(X^{*}))|X,D=1]$ and $\Pr[D=1|X]$.

##  Appendix C Additional results from transportability analyses using NLST
and NHANES data

Table S2 shows characteristics of the variables used in the NLST to NHANES
analysis presented in Section 5 in the main paper stratified by datasource.

Table S2: Summary of the variables used in the transportability analysis of
lung cancer prediction models in Section 5 of the main paper. Continuous
variables are summarized by their weighted mean (weighted standard deviation)
and categorical variables are summarized by weighted percentage in each
category.

[htbp]  Variable NLST NHANES Age 61.4 (5.0) 62.6(5.4) BMI 27.9 (5.1) 29.1
(6.5) Race or ethnic group Black 4.4% 8.2% White 89.9% 81.6% Hispanic 1.8%
6.2% Other 3.8% 4.0% Education level less than high school 6.3% 36.2% high
school graduate 38.3% 30.3% AA degree/some college 23.7% 21.6% college
graduate 31.7% 11.9% Personal history of cancer 4.1% 22.3% Smoking status
Current 48.4% 62.4% Former 51.6% 37.6% Smoking intensity 11footnotemark: 1
28.5 (11.5) 27.2(12.1) Duration of smoking (year) 39.8 (7.3) 42.1 (6.7)
Smoking quit time (year) 22footnotemark: 2 7.3 (4.8) 6.4 (3.7)

  * 1

Smoking intensity (the average number of cigarettes smoked per day) has a
nonlinear association with lung cancer and was transformed by dividing by 10,
exponentiating by the power -1 in the prediction model for lung cancer.

  * 2

Smoking quit time in former smokers (years).

##  Appendix D Evaluation of a prediction model for a positive lung screen

In this section, we used only data from the National Lung Screening Trial
(NLST) to evaluate the performance of the different estimators for the target
population risk $\psi$. Here we focused on the binary outcome of whether the
baseline $T_{0}$ low-dose CT screen was positive (the rate of positive screens
was $27.4\%$) and only used data from the CT arm of the NLST trial.

To compare the different estimators of the Brier risk in the target
population, we simulated whether an observation was from the source population
using two settings. In the first setting, the probability of being from the
source population was $0.5$ for all covariate patterns, which implies that the
covariate distribution was the same in the source population and the target
population (and hence the naïve empirical estimator is expected to be
unbiased). In the second setting, the covariate distributions differed between
the two populations. The probability of being from the source population was
simulated using a main effects logistic regression model with coefficients
equal to $0.05$ for all covariates that had a positive association with the
outcome (determined by fitting a main effects logistic regression model) and
$-0.05$ for all covariates that had a negative association with the outcome
(introducing systematic differences between the sample from the source
population and the sample from the target population). For both settings, we
simulated 1000 different splits into source and target population data.

For both settings, the prediction model, $g(X^{*})$, we assessed the
performance of was a main effects logistic regression model that included age
at randomization, height, weight, pack years of smoking, average number of
cigarettes smoked per day, age started smoking, total years of smoking, race,
gender, education status, marital status, whether has smoked cigars, whether
has smoked pipes, whether lives with a smoker, whether works with a smoker,
whether has worked with asbestos for at least one year, whether has worked on
a farm for at least one year, prior diagnosis of emphysema, prior diagnosis of
hypertension, prior diagnosis of pneumonia, whether siblings have had lung
cancer, and whether parents have had lung cancer. And those were also the
variables we adjusted for in the transportability analysis (i.e., $X^{*}=X)$.
The prediction model, $g(X^{*})$, was fit using the source data training set,
which was a random sample of 2/3 of the source data. Model performance was
evaluated using the Brier loss estimated using the combination of the test set
from the source data and all of the target data. For all the estimators,
standard deviation was estimated using the non-parametric bootstrap with 1000
bootstrap samples.

Both estimators $\widehat{p}(X)$ and $\widehat{h}(X)$ were from a linear main
effects logistic regression model. As outcome information was available in the
target data, we can get an unbiased estimator for
$\operatorname{\textnormal{\mbox{E}}}[(Y-g(X^{*}))^{2}|D=0]$ by using only
target population data. This target population only estimator will be treated
as an “oracle” estimator and used for comparisons to evaluate the bias
(compared to the “oracle” estimator) of the four estimators of $\psi$.

Table 3(b) shows results for both settings averaged across $1000$ different
splits into target and source data. The estimates from the oracle estimator
that used outcome data from the simulated target population was 0.197
(standard error: 0.00188) for the setting when the covariate distributions are
the same in the target population and the source population and was 0.179 (SD:
0.00654) when the covariate distributions differ between the two populations.

For the setting where the covariate distribution was the same in the source
and the target data, all four estimates were unbiased and had almost identical
estimated standard deviation. For the setting where the covariate
distributions differed, the estimates from the naïve empirical estimator were
biased (11% relative bias) while the estimates from the other three estimators
were approximately unbiased. Among the unbiased estimators, the conditional
loss estimator had the smallest estimated standard deviation, followed by the
doubly robust estimator, and the weighting estimator had the largest estimated
standard deviation.

Table S3: Comparison of the naïve empirical source estimator, weighting
estimator, conditional loss estimator, and the doubly robust estimator using
data from the National Lung Cancer Screening Trial. Bias is the difference
between the average estimates from the estimator and the oracle target
population estimator that used outcome information from the target population.
Standard deviation is the average estimated standard deviation of the
estimator (estimated using the non-parametric bootstrap with 1000 bootstrap
samples). The results presented are averaged across 1000 simulated splits into
a source and target data. Table (a) shows results from the simulation setting
where covariate distributions in the source and target data were the same and
Table (b) shows results when the two covariate distributions differed. To help
with numerical comparisons, we multiplied the simulation estimates of both the
bias and average standard deviation by $10^{3}$.

(a) Covariate distributions are the same

Bias Standard deviation Naïve empirical 0.0642 3.19 Weighting 0.0514 3.20
Conditional loss 0.0540 3.19 Doubly robust 0.0525 3.20

(b) Covariate distributions differ

Bias Standard deviation Naïve empirical 18.8 2.29 Weighting 0.754 23.2
Conditional loss -0.498 57.6 Doubly robust 1.43 21.2