###### Proof.

.: 5pt

######  Systematically Missing Data in Causally Interpretable Meta-Analysis

###### Abstract

Causally interpretable meta-analysis combines information from a collection of
randomized controlled trials to estimate treatment effects in a target
population in which experimentation may not be possible but covariate
information can be collected from a simple random sample. In such analyses, a
key practical challenge is systematically missing data when some baseline
covariates are not collected in all trials. Here, we provide identification
results for potential (counterfactual) outcome means and average treatment
effects in the target population when covariate data are systematically
missing from some of the trials in the meta-analysis. We propose three
estimators for the average treatment effect in the target population, examine
their asymptotic properties, and show that they have good finite-sample
performance in simulation studies. We use the estimators to analyze data from
two large lung cancer screening trials and target population data from the
National Health and Nutrition Examination Survey (NHANES). To accommodate the
complex survey design of the NHANES, we modify the methods to incorporate
survey sampling weights and allow for clustering.

Keywords: causally interpretable meta-analysis, domain adaptation,
systematically missing data, transportability, generalizability, covariate
shift, multi-source inference

##  1 Introduction

When multiple randomized trials compare the same treatments, it is natural to
want to learn about treatment effects by synthesizing evidence across trials.
Meta-analysis is an umbrella term for quantitative methods for evidence
synthesis [1, 2]. When each trial included in a meta-analysis recruits
participants from a different, typically ill-defined, underlying population
and treatment effects are heterogeneous across populations, the “summary”
estimates produced by standard meta-analysis methods do not have a clear
causal interpretation because they cannot be interpreted as treatment effects
in a well-defined target population [3]. Building on methods for extending
(generalizing or transporting [4]) treatment effects from a single clinical
trial to a target population of substantive interest [5, 6, 7, 8, 9], we
recently proposed methods for “causally interpretable meta-analysis” that
combines information from multiple trials to estimate treatment effects in a
well-defined target population [10, 3].

These methods assume that a common set of covariates adequate to render the
trials and the target population exchangeable is available from all the trials
and the sample of the target population. In attempts to apply the methods in
practice, we have found that this assumption is often not true [11] and
instead one or more covariates are systematically missing, in the sense that
they are available from some trials but not others. This situation is
different from the usual within-trial missingness, where values are missing
from some observations in a given trial but available from other observations
in the same trial: with systematic missingness, for one or more covariates,
values are missing from all observations in a given trial.

There is a large literature on methods for dealing with within-trial missing
data [12, 13, 14]. Popular methods for handling within-trial missing data
assume that there is a positive probability of observing the full data given
any values of the observed data. If a covariate is not measured in a specific
data source, the probability of observing the full data conditional on an
observation being from that data source is zero, violating the assumptions
needed for within-trial missing data methods. Thus, addressing systematically
missing data requires different identifiability conditions and estimation
procedures than those used for within-trial missing data. Previous work on
systematically missing data has focused on conventional individual participant
data meta-analysis methods and has only considered multiple imputation-based
approaches for addressing missingness [15, 16, 17, 18, 19]. The approaches
allow all covariates collected in each trial to be used for analysis and
therefore can account for variables that are collected differently across
trials (e.g., different trials use different survey instruments or measurement
devices). This has advantages over complete case analysis that either does not
use trials with missing covariates (leading to loss in efficiency) or does not
use covariates that are missing in some trials (leading to potential bias).
Despite these advantages, however, existing methods have not been given a
clear causal interpretation, are often only supported by intuitive arguments
(e.g., without formal identification analysis), and their theoretical
properties remain largely unknown. In fact, the vast majority of the
imputation methods rely on chained equations where the final imputation model
may not even correspond to a true joint density [20, 21].

Here, we propose methods for handling systematically missing data when
extending inferences about potential (counterfactual) outcome means and
average treatment effects from multiple trials to a target population that
lacks outcome or treatment information. We provide identification results,
propose g-formula, weighting, and augmented weighting (doubly robust)
estimators for these causal estimands, and study large-sample behavior of the
estimators. We also examine the finite-sample behavior of the estimators using
simulations and implement them to analyze data from two trials of lung cancer
screening to learn about the effects of the treatments on a nationally
representative target population. The target population data are obtained from
the National Health and Nutrition Examination Survey (NHANES) and we show how
to modify our estimators to account for the NHANES complex survey design by
incorporating the survey sampling weights and multi-stage clustering.

##  2 Data structure and causal estimands

We have data from a collection of studies $\mathcal{S}$, indexed by
$s=1,\ldots,L$, evaluating the effect of the same treatments on the same
outcome. Although we will refer to these data sources as “trials” in the
remainder of the paper, we note here that the collection $\mathcal{S}$ can
consist of observational studies or a mixture of trials and observational
studies (provided the conditions stated below can be plausibly assumed to hold
for them). Let $Y$ denote a univariate outcome (continuous, binary, or count)
measured at the end of followup, $A$ the random variable for treatment
assignment taking values in a finite set $\mathcal{A}$, $X$ the vector of
baseline (pre-randomization and pre-treatment) covariates, and $S$ a random
variable indicating which trial the data comes from. In the absence of any
missing data, the data from trial $s\in\mathcal{S}$ is assumed to be
$\\{(X_{i},A_{i},Y_{i},S_{i}=s):i=1,\ldots,n_{s}\\}$, where $n_{s}$ is the
total number of observations in trial $s$.

We also collect baseline covariate information from a separately obtained
random sample of $n_{0}$ individuals from the target population of substantive
interest: $\\{X_{i}:i=1,\ldots,n_{0}\\}$. We shall use $S=0$ to indicate
observations from the target population. We do not require treatment or
outcome information from the sample of the target population – this is a
possible strength of our methods because they can be used when treatment and
outcome data from the target population are not collected (e.g., due to cost
or the need for specialized ascertainment procedures), are of insufficient
quality (e.g., due to gross measurement error), or subject to confounding by
unmeasured variables (e.g., because treatment is not randomly assigned in the
target population). Let $n=\sum_{j=0}^{L}n_{j}$ be the total sample size for
the composite dataset formed by appending data from the collection of trials
and the sample of the target population. In the absence of missing data, the
data available on participant $i\in\\{1,\ldots,n\\}$ would be a realization of

$\mathcal{F}_{i}=\begin{cases}(X_{i},A_{i},Y_{i},S_{i}),&\text{if
$S_{i}\in\mathcal{S}$},\\\ (X_{i},S_{i}=0),&\text{if $S_{i}=0$}.\end{cases}$
(1)

Hereafter, we refer to $\mathcal{F}=\\{\mathcal{F}_{i}:i=1,\ldots,n\\}$ as the
full data. Let $Y_{i}^{a}$ denote the potential outcome for participant $i$
under intervention to set treatment to $a\in\mathcal{A}$ [22, 23]. We focus on
estimating the potential outcome means in the target population
$\operatorname{\textnormal{\mbox{E}}}[Y^{a}|S=0]$, for $a\in\mathcal{A}$. Many
causal estimands of interest are functions of potential outcome means; for
example the average treatment effect comparing two treatments $a$ and
$a^{\prime}$ in the target population is equal to the difference of the
corresponding potential outcome means:
$\operatorname{\textnormal{\mbox{E}}}[Y^{a}-Y^{a^{\prime}}|S=0]=\operatorname{\textnormal{\mbox{E}}}[Y^{a}|S=0]-\operatorname{\textnormal{\mbox{E}}}[Y^{a^{\prime}}|S=0]$.

Suppose now that some covariates are not collected in some of the trials
(i.e., one or more covariates are systematically missing from some trials in
the collection $\mathcal{S}$) but all covariates are collected in the target
population. Let $K$ denote the number of different missing data patterns of
systematically missing covariates; $X^{(k)}$ the observed covariate vector
under missingness pattern $k\in\\{1,\ldots,K\\}$; $X^{(-k)}$ the a vector of
the the components of $X$ that are missing under missingness pattern
$k\in\\{1,\ldots,K\\}$; and $\mathcal{S}^{k}$ be the set of trials with
missingness pattern $k$. Note that we allow the systematic missingness to be
potentially non-monotone across trials. The observed data on participant
$i\in\\{1,\ldots,n\\}$ are a realization of

$\mathcal{O}_{i}=\begin{cases}(X_{i}^{(k)},A_{i},Y_{i},S_{i}),&\text{if
$S_{i}\in\mathcal{S}^{k}$},\\\ (X_{i},S_{i}=0),&\text{if
$S_{i}=0$}.\end{cases}$ (2)

Define $n_{k}^{*}=\sum_{i=1}^{n}I(S_{i}\in\mathcal{S}^{k})$ as the number of
observations with missingness pattern $\mathcal{S}^{k}$ and assume that
$\frac{n_{k}^{*}}{n}\longrightarrow p_{k}>0$ as $n\longrightarrow\infty$. We
use the notation $\Pr_{k}$ and $\operatorname{\textnormal{\mbox{E}}}_{k}$ to
denote probabilities and expectations, respectively, that are conditional on
the data sources in $\mathcal{S}^{k}\cup\\{0\\}$. Throughout, we use
$f(\cdot)$ to generically denote densities.

##  3 Identification analysis

###  3.1 Identifiability conditions

We will show that the potential outcome means in the target population are
identifiable using the observable data under the following identifiability
conditions:

  1. A1.

Consistency: If $A_{i}=a$, then $Y_{i}^{a}=Y_{i}$ for every individual $i$ in
any of the trials or the target population.

  2. A2.

Within-trial exchangeability over treatment for trials with missingness
pattern $k$: $Y^{a}\mathbin{\mathchoice{\hbox
to0.0pt{\hbox{\set@color$\displaystyle\perp$}\hss}\kern 3.46875pt{}\kern
3.46875pt\hbox{\set@color$\displaystyle\perp$}}{\hbox
to0.0pt{\hbox{\set@color$\textstyle\perp$}\hss}\kern 3.46875pt{}\kern
3.46875pt\hbox{\set@color$\textstyle\perp$}}{\hbox
to0.0pt{\hbox{\set@color$\scriptstyle\perp$}\hss}\kern 2.36812pt{}\kern
2.36812pt\hbox{\set@color$\scriptstyle\perp$}}{\hbox
to0.0pt{\hbox{\set@color$\scriptscriptstyle\perp$}\hss}\kern 1.63437pt{}\kern
1.63437pt\hbox{\set@color$\scriptscriptstyle\perp$}}}A|(X,S,S\in\mathcal{S}^{k})$
for each $k\in\\{1,\ldots,K\\}$ and $a\in\mathcal{A}$.

  3. A3.

Positivity of treatment assignment for trials with missingness pattern $k$:
For every treatment $a\in\mathcal{A}$ and all missingness patterns
$k\in\\{1,\ldots,K\\}$, $\Pr_{k}[A=a|X^{(k)}=x^{(k)},S\in\mathcal{S}^{k}]>0$
for every $x^{(k)}$ such that $f(x^{(k)},S\in\mathcal{S}^{k})>0$.

  4. A4.

Exchangeability over data source (transportability):
$Y^{a}\mathbin{\mathchoice{\hbox
to0.0pt{\hbox{\set@color$\displaystyle\perp$}\hss}\kern 3.46875pt{}\kern
3.46875pt\hbox{\set@color$\displaystyle\perp$}}{\hbox
to0.0pt{\hbox{\set@color$\textstyle\perp$}\hss}\kern 3.46875pt{}\kern
3.46875pt\hbox{\set@color$\textstyle\perp$}}{\hbox
to0.0pt{\hbox{\set@color$\scriptstyle\perp$}\hss}\kern 2.36812pt{}\kern
2.36812pt\hbox{\set@color$\scriptstyle\perp$}}{\hbox
to0.0pt{\hbox{\set@color$\scriptscriptstyle\perp$}\hss}\kern 1.63437pt{}\kern
1.63437pt\hbox{\set@color$\scriptscriptstyle\perp$}}}S|X$ for each
$a\in\mathcal{A}$.

  5. A5.

Positivity of participation for trials for a given missingness pattern:
$\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}=x^{(k)}]>0$ for each
$k\in\\{1,\ldots,K\\}$ and for every $x^{(k)}$ that has a positive density in
the target population, $f(x^{(k)},S=0)>0$.

  6. A6.

$Y\mathbin{\mathchoice{\hbox
to0.0pt{\hbox{\set@color$\displaystyle\perp$}\hss}\kern 3.46875pt{}\kern
3.46875pt\hbox{\set@color$\displaystyle\perp$}}{\hbox
to0.0pt{\hbox{\set@color$\textstyle\perp$}\hss}\kern 3.46875pt{}\kern
3.46875pt\hbox{\set@color$\textstyle\perp$}}{\hbox
to0.0pt{\hbox{\set@color$\scriptstyle\perp$}\hss}\kern 2.36812pt{}\kern
2.36812pt\hbox{\set@color$\scriptstyle\perp$}}{\hbox
to0.0pt{\hbox{\set@color$\scriptscriptstyle\perp$}\hss}\kern 1.63437pt{}\kern
1.63437pt\hbox{\set@color$\scriptscriptstyle\perp$}}}X^{(-k)}|X^{(k)},S\in\mathcal{S}^{k},A=a$
for every treatment $a\in\mathcal{A}$ and each missingness pattern
$k\in\\{1,\ldots,K\\}$.

Condition A1 implies that there is no direct effect of participation in some
specific trial on the outcome [24]. Conditions A2 and A3 are expected to hold
when the treatment is randomized and are on the collection of trials with
missingness pattern $k$. Conditions A4 and A5 are used for causally
interpretable meta-analysis in the absence of systematically missing data [10]
and condition A6 ensures that for each missingness pattern the observable data
are sufficient to identify the potential outcome means in the target
population. We note that positivity of trial participation (condition A5) is
on the pooled set from all trials with missingness pattern $k$. This implies
that we can draw inference about a target population that has a broader
covariate distribution than any individual trial. For example, assume a
missingness pattern has two trials and the only covariate needed for
transportability to the target population is age. If the first trial enrolls
people aged $50-75$ and the second trial enrolls people aged $40-70$, then the
positivity condition is satisfied for a target population that has an age
range of $40-75$. Similarly, the positivity of treatment assignment condition
(A3) is weaker than assuming positivity of treatment assignment within each
trial as it only requires at least one trial in $\mathcal{S}^{k}$ to have a
positive probability of assigning participants with missingness pattern $k$ to
treatment $a$. The conditions imply that for a fixed treatment
$a\in\mathcal{A}$, the expectations
$\operatorname{\textnormal{\mbox{E}}}[Y|X^{(k)},S=s,A=a]$ are equal for all
trials $s\in\mathcal{S}^{k}$ [10]. Because the expectations
$\operatorname{\textnormal{\mbox{E}}}[Y|X^{(k)},S=s,A=a]$ are identified from
the observable data for each $s\in\mathcal{S}^{k}$, the equality of these
expectations is testable. Throughout, to simplify exposition, we assume
perfect adherence to the the assigned treatment strategy, no losses to
followup, and no missing data. In practical applications, when these
complications arise, our methods can be naturally combined with well-
established approaches for addressing them.

###  3.2 Identification

In Supplementary Web Appendix A.1, we prove the following identification
result:

######  Theorem 1.

Under Assumptions A1 through A6, for each missing data pattern
$k\in\\{1,\ldots,K\\}$, the potential outcome mean in the target population,
$\operatorname{\textnormal{\mbox{E}}}[Y^{a}|S=0]$, can be written as the
following observed data functional:

$\psi_{k}(a)=\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},A=a,S\in\mathcal{S}^{k}]|S=0],$
(3)

or, equivalently, using the weighting representation

$\psi_{k}(a)=\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\left[\frac{\Pr_{k}[S=0|X^{(k)}]I(S\in\mathcal{S}^{k},A=a)}{\Pr_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]}Y\right].$
(4)

From results (3) and (4), we have the following corollary.

######  Corollary 1.

For all weight vectors $\\{w_{k}:k=1,\ldots,K\\}$ satisfying
$\sum_{k=1}^{K}w_{k}=1$, the potential outcome mean in the target population
can be identified by

$\sum_{k=1}^{K}w_{k}\psi_{k}(a)=\sum_{k=1}^{K}w_{k}\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},A=a,S\in\mathcal{S}^{k}]|S=0],$

or, equivalently, by the weighting expression

$\sum_{k=1}^{K}w_{k}\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\left[\frac{\Pr_{k}[S=0|X^{(k)}]I(S\in\mathcal{S}^{k},A=a)}{\Pr_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]}Y\right].$

We discuss the choice of the weights $\\{w_{1},\ldots,w_{k}\\}$ in Section 4.

The sampling model for causally interpretable meta-analysis assumes a
stratified sampling setup, where trial observations are treated as-if randomly
sampled from a (typically unspecified) super-population stratified by $S$ with
sampling fractions that are constant within a trial, can vary between trials,
and are unknown to the investigators [3]. Furthermore, we assume that the
target population data are randomly sampled from the target population,
separately from the trial samples and with unknown sampling fraction. When the
sampling fractions vary and are unknown, the probabilities
$\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]$ are not identifiable as the sample
proportion of observations from $\mathcal{S}^{k}$ in the data does not in
general reflect the probability of being from $\mathcal{S}^{k}$ in the super-
population. Nevertheless, as we show in Supplementary Web Appendix B, the
potential outcome mean is still identifiable. More specifically, we show that
if $D$ is an indicator that is equal to one when an individual contributes
data to the analyses and zero otherwise, the potential outcome mean can be
written as

$\psi_{k}(a)=\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},A=a,S\in\mathcal{S}^{k},D=1]|S=0,D=1],$
(5)

or using the weighting representation

$\psi_{k}(a)=\frac{1}{\Pr_{k}[S=0|D=1]}\operatorname{\textnormal{\mbox{E}}}_{k}\left[\frac{\Pr_{k}[S=0|X^{(k)},D=1]I(S\in\mathcal{S}^{k},A=a)}{\Pr_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k},D=1]\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)},D=1]}Y\Bigg{|}D=1\right].$
(6)

Because all quantities in expressions (5) and (6) condition on data
availability ($D=1$), the potential outcome means are still identifiable under
the sampling model described above. It follows that the potential outcome
means and treatments effects in the target population are identifiable under
this sampling model and the expectations in Theorem 1 can be interpreted as
integrals with respect to the densities under the same sampling model.

##  4 Estimation and inference

###  4.1 Estimation

Up to now, we have treated the weights $w_{k}$ as fixed. But they can be
estimated using the data; we use $\widehat{w}_{k}$ to denote the potentially
data-dependent weights. We denote the asymptotic limit of $\widehat{w}_{k}$ as
$w^{*}_{k}>0$ for $k=1,\ldots,K$ and we assume that, for each $k$,
$\widehat{w}_{k}$ is $\sqrt{n}$-convergent, and
$\sum_{k=1}^{K}\widehat{w}_{k}=1$. By using plug-in estimators for the unknown
quantities in expressions (3) and (4) we obtain two estimators for the
potential outcome mean in the target population. The g-formula [25] (outcome
model-based) estimator is

$\widehat{\psi}_{\text{\tiny
g}}(a)=\sum_{k=1}^{K}\widehat{w}_{k}\frac{1}{n_{0}}\sum_{i=1}^{n}I(S_{i}=0)\widehat{g}_{a,k}(X_{i}^{(k)}).$
(7)

Here, $\widehat{g}_{a,k}(X^{(k)})$ is an estimator for
$\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},A=a,S\in\mathcal{S}^{k}]$.
The g-formula estimator $\widehat{\psi}_{{\tiny g}}(a)$ is consistent if all
outcome models are correctly specified, that is, if
$\widehat{g}_{a,k}(X^{(k)})\overset{P}{\longrightarrow}\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},A=a,S\in\mathcal{S}^{k}]$
for each $k=1,\ldots,K$.

The weighting estimator for the potential outcome mean in the target
population is

$\widehat{\psi}_{\text{\tiny
w}}(a)=\sum_{k=1}^{K}\widehat{w}_{k}\frac{1}{n_{0}}\sum_{i=1}^{n}\widehat{o}_{a}(X^{(k)}_{i},A_{i},S_{i})Y_{i},$
(8)

where, for each missingness pattern $k\in\\{1,\ldots,K\\}$, we define the
weights $\widehat{o}_{a}(X^{(k)},A,S)$ as

$\widehat{o}_{a}(X^{(k)},A,S)=\frac{(1-\widehat{p}_{k}(X^{(k)}))I(S\in\mathcal{S}^{k},A=a)}{\widehat{e}_{a,k}(X^{(k)})\widehat{p}_{k}(X^{(k)})}.$

Here, $\widehat{p}_{k}(X^{(k)})$ is an estimator for
$\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]$ and $\widehat{e}_{a,k}(X^{(k)})$ is an
estimator for $\Pr_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]$. The weighting
estimator is consistent if both the models for $\widehat{e}_{a,k}(X^{(k)})$
and $\widehat{p}_{k}(X^{(k)})$ are correctly specified, that is, if
$\widehat{e}_{a,k}(X^{(k)})\overset{P}{\longrightarrow}\Pr_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]$
and
$\widehat{p}_{k}(X^{(k)})\overset{P}{\longrightarrow}\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]$,
for each $k=1,\ldots,K$.

The weighting estimator is not guaranteed to take values in the support of
$Y$, but a “normalized” version [26, 27] of the weighting estimator

$\widehat{\psi}_{\text{\tiny
w,norm}}(a)=\sum_{k=1}^{K}\widehat{w}_{k}\frac{1}{\sum_{i=1}^{n}\widehat{o}_{a}(X^{(k)}_{i},A_{i},S_{i})}\sum_{i=1}^{n}\widehat{o}_{a}(X^{(k)}_{i},A_{i},S_{i})Y_{i}$

is guaranteed to take values in the support of $Y$.

In Supplementary Web Appendix A.2, we derive the first-order influence
function of $\psi(a)$ and using plug-in estimators into the unknown quantities
results in the “doubly robust” [28] estimator

$\widehat{\psi}_{\text{\tiny
DR}}(a)=\sum_{k=1}^{K}\widehat{w}_{k}\frac{1}{n_{0}}\sum_{i=1}^{n}\left(\widehat{o}_{a}(X^{(k)}_{i},A_{i},S_{i})(Y_{i}-\widehat{g}_{a,k}(X^{(k)}_{i}))+I(S_{i}=0)\widehat{g}_{a,k}(X^{(k)}_{i})\right).$
(9)

A normalized version of the doubly robust estimator is given by

$\widehat{\psi}_{\text{\tiny
DR,norm}}(a)=\sum_{k=1}^{K}\widehat{w}_{k}\left(\frac{\sum_{i=1}^{n}\widehat{o}_{a}(X^{(k)}_{i},A_{i},S_{i})(Y_{i}-\widehat{g}_{a,k}(X^{(k)}_{i}))}{\sum_{i=1}^{n}\widehat{o}_{a}(X^{(k)}_{i},A_{i},S_{i})}+\frac{1}{n_{0}}\sum_{i=1}^{n}I(S_{i}=0)\widehat{g}_{a,k}(X^{(k)}_{i})\right).$

The normalized version of the doubly robust estimator is not guaranteed to
take values in the support of $Y$, but will often have better finite sample
performance when some weights are extreme [29, 30]. The g-formula and the
weighting estimators are special cases of the doubly robust estimator obtained
by setting $\widehat{p}_{k}(X^{(k)})=1$ and $\widehat{g}_{a,k}(X^{(k)})=0$,
respectively. All three estimators $\widehat{\psi}_{{\tiny
g}}(a),\widehat{\psi}_{\text{\tiny w}}(a),$ and $\widehat{\psi}_{\text{\tiny
DR}}(a)$ can be calculated using the observed data as they do not require
outcome or treatment information in the target population and for trials in
the set of trials with missingness pattern $k$, that is $\mathcal{S}^{k}$, the
estimators only depend on the observed covariate pattern $X^{(k)}$.

###  4.2 Large-sample behavior of the doubly robust estimator

We define some notation that will be useful in deriving the large-sample
behavior of $\widehat{\psi}_{\text{\tiny DR}}(a)$. For each
$k\in\\{1,\ldots,K\\}$ and general functions
$g_{a,k}^{\prime}(X^{(k)}),e_{a,k}^{\prime}(X^{(k)}),\gamma_{k}^{\prime},$ and
$p^{\prime}_{k}(X^{(k)})$, define

$\displaystyle L(g_{a,k}^{\prime}(X^{(k)}),$ $\displaystyle
e_{a,k}^{\prime}(X^{(k)}),\gamma_{k}^{\prime},p^{\prime}_{k}(X^{(k)}))$
$\displaystyle=\gamma_{k}^{\prime}\left(\frac{(1-p^{\prime}_{k}(X^{(k)}))I(S\in\mathcal{S}^{k},A=a)}{e_{a,k}^{\prime}(X^{(k)})p^{\prime}_{k}(X^{(k)})}(Y-g_{a,k}^{\prime}(X^{(k)}))+I(S=0)g_{a,k}^{\prime}(X^{(k)})\right).$

Define $\widehat{\gamma}_{k}=\frac{n_{0}+n_{k}}{n_{0}}$, which is a non-
parametric estimator for $\gamma_{k}=\Pr_{k}[S=0]^{-1}$. For a random variable
$W$ define

$\mathbb{P}_{k}(W)=\frac{1}{n_{k}+n_{0}}\sum_{i:S_{i}\in\\{0\\}\cup\mathcal{S}^{k}}W_{i}$

and
$\mathbb{G}_{k}(W)=\sqrt{n}(\mathbb{P}_{k}(W)-\operatorname{\textnormal{\mbox{E}}}_{k}[W])$.
Using these definitions, we can rewrite the estimator as

$\widehat{\psi}_{\text{\tiny
DR}}(a)=\sum_{k=1}^{K}\widehat{w}_{k}\mathbb{P}_{k}\big{(}L(\widehat{g}_{a,k}(X^{(k)}),\widehat{e}_{a,k}(X^{(k)}),\widehat{\gamma}_{k},\widehat{p}_{k}(X^{(k)}))\big{)}.$

Let $g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}_{k}(X^{(k)})$
be the asymptotic limits of
$\widehat{g}_{a,k}(X^{(k)}),\widehat{e}_{a,k}(X^{(k)}),\widehat{\gamma}_{k},\widehat{p}_{k}(X^{(k)})$,
respectively. We assume that all the limits
$g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}_{k}(X^{(k)})$
exist and

$\mathbb{P}_{k}\big{(}L(\widehat{g}_{a,k}(X^{(k)}),\widehat{e}_{a,k}(X^{(k)}),\widehat{\gamma}_{k},\widehat{p}_{k}(X^{(k)}))\big{)}$

is uniformly bounded. To derive the asymptotic distribution of
$\widehat{\psi}_{\text{\tiny DR}}(a)$ we make the following assumptions:

  1. B1.

For each $k\in\\{1,\ldots,K\\}$, at least one of the two following conditions
hold

$\widehat{g}_{a,k}(X^{(k)})\overset{P}{\longrightarrow}\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},S\in\mathcal{S}^{k},A=a]$

or

$\displaystyle\widehat{p}_{k}(X^{(k)})$
$\displaystyle\overset{P}{\longrightarrow}{\Pr}_{k}[S\in\mathcal{S}^{k}|X^{(k)}]$
$\displaystyle\widehat{e}_{a,k}(X^{(k)})$
$\displaystyle\overset{P}{\longrightarrow}{\Pr}_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}].$

  2. B2.

For each $k\in\\{1,\ldots,K\\}$, the sequences
$L\big{(}\widehat{g}_{a,k}(X^{(k)}),\widehat{e}_{a,k}(X^{(k)}),\widehat{\gamma}_{k},\widehat{p}_{k}(X^{(k)})\big{)}$
and
$L\big{(}g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}_{k}(X^{(k)})\big{)}$
are Donsker [31].

  3. B3.

$||L\big{(}\widehat{g}_{a,k}(X^{(k)}),\widehat{e}_{a,k}(X^{(k)}),\widehat{\gamma}_{k},\widehat{p}_{k}(X^{(k)})\big{)}-L\big{(}g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}_{k}(X^{(k)})\big{)}||_{2}\overset{P}{\longrightarrow}0.$

  4. B4.

$\operatorname{\textnormal{\mbox{E}}}\big{[}L\big{(}g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}_{k}(X^{(k)})\big{)}^{2}\big{]}<\infty$.

Assumption B1 reflects the model doubly robustness property as the estimator
requires for each $k\in\\{1,\ldots,K\\}$ that at least one of the outcome
model $\widehat{g}_{a,k}(X^{(k)})$ or both the model for trial participation
$\widehat{p}_{k}(X^{(k)})$ and the treatment model
$\widehat{e}_{a,k}(X^{(k)})$ are correctly specified. Note that different
missingness patterns can satisfy different parts of assumption B1. In other
words, the assumption can be satisfied if for some patterns we are only able
to correctly specify the models for the outcome and for others we are only
able to correctly specify the models for trial participation and treatment,
provided that at least one of these groups of models is satisfied for each
pattern. Assumption B2 follows from Donsker preservation theorems if all the
individual estimators
$\widehat{g}_{a,k}(X^{(k)}),\widehat{e}_{a,k}(X^{(k)}),\widehat{p}_{k}(X^{(k)})$
and their corresponding limits are Donsker, the models
$\widehat{e}_{a,k}(X^{(k)})$ and $\widehat{p}_{k}(X^{(k)})$ are uniformly
bounded away from zero, and $\widehat{g}_{a,k}(X^{(k)})$ is uniformly bounded.
The Donsker assumption restricts the class of models, but includes many
commonly used classes of models that are more flexible than standard
parametric models [32].

######  Theorem 2.

If assumptions B1 through B4 hold, then:

  1. 1.

The doubly robust estimator is consistent,
$\widehat{\psi}_{\text{\tiny\emph{DR}}}(a)\overset{P}{\longrightarrow}\psi(a).$

  2. 2.

The doubly robust estimator has the asymptotic representation

$\sqrt{n}\big{(}\widehat{\psi}_{\text{\tiny\emph{DR}}}(a)-\psi(a)\big{)}=\sum_{k=1}^{K}w^{*}_{k}\mathbb{G}_{k}\bigg{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}_{k}(X^{(k)}))\bigg{)}+REM+o_{p}(1),$
(10)

where

$\displaystyle REM\leq$
$\displaystyle\sqrt{n}O_{P}\Bigg{(}\sum_{k=1}^{K}(\widehat{w}_{k}-w^{*}_{k})+\sum_{k=1}^{K}\Big{|}\Big{|}\widehat{g}_{a,k}(X^{(k)})-\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},S\in\mathcal{S}^{k},A=a]\Big{|}\Big{|}_{2}$
$\displaystyle\times\Big{(}\Big{|}\Big{|}\widehat{p}_{k}(X^{(k)})-{\Pr}_{k}[S\in\mathcal{S}^{k}|X^{(k)}]\Big{|}\Big{|}_{2}+\Big{|}\Big{|}\widehat{e}_{a,k}(X^{(k)})-{\Pr}_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]\Big{|}\Big{|}_{2}\Big{)}\Bigg{)}.$

The second part of Theorem 2 gives the rate of convergence result

$\displaystyle||\widehat{\psi}_{\text{\tiny DR}}(a)$
$\displaystyle-\psi(a)||_{2}\leq
O_{P}\Bigg{(}\frac{1}{\sqrt{n}}+\sum_{k=1}^{K}\Big{|}\Big{|}\widehat{g}_{a,k}(X^{(k)})-\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},S\in\mathcal{S}^{k},A=a]\Big{|}\Big{|}_{2}$
$\displaystyle\times\Big{(}\Big{|}\Big{|}\widehat{p}_{k}(X^{(k)})-{\Pr}_{k}[S\in\mathcal{S}^{k}|X^{(k)}]\Big{|}\Big{|}_{2}+\Big{|}\Big{|}\widehat{e}_{a,k}(X^{(k)})-{\Pr}_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]\Big{|}\Big{|}_{2}\Big{)}\Bigg{)}.$

This shows that the rate of convergence of the doubly robust estimator depends
on how fast the rate of convergence of the nuisance estimators
$\widehat{g}_{a,k}(X^{(k)}),\widehat{p}_{k}(X^{(k)}),$ and
$\widehat{e}_{a,k}(X^{(k)})$ are. For example, if for each
$k\in\\{1,\ldots,K\\}$ the combined rate of convergence of
$\widehat{g}_{a,k}(X^{(k)})$ and
$(\widehat{p}_{k}(X^{(k)}),\widehat{e}_{a,k}(X^{(k)}))$ is $\sqrt{n}$, then
the rate of convergence for $\widehat{\psi}_{\text{\tiny DR}}(a)$ is
$\sqrt{n}$. This allows the nuisance parameters to be estimated using
procedures that converge at rates slower than $\sqrt{n}$, as long as the
combined rate of convergence is $\sqrt{n}$ (rate double robustness [33]). This
is not the case for the outcome or weighting estimators that inherit the rate
of convergence of their respective nuisance parameters,
$\widehat{g}_{a,k}(X^{(k)})$, or
$(\widehat{p}_{k}(X^{(k)}),\widehat{e}_{a,k}(X^{(k)}))$.

Although the Donsker assumption on the estimators for the nuisance parameters
(Assumption B2) is more flexible than assuming a parametric model, it can be
too restrictive for more data-adaptive estimators. In that case, sample
splitting [34, 35, 36] can be used to avoid any Donsker assumptions. The
sample split doubly robust estimator is calculated by:

  1. 1.

Split the data into two mutually exclusive and exhaustive sets. Let
$\mathcal{B}_{1}$ and $\mathcal{B}_{2}$ be indexes from each set.

  2. 2.

For each missingness pattern, use observations with indexes in
$\mathcal{B}_{2}$ to calculate the estimators
$\widehat{g}_{a,k}^{(2)}(X^{(k)}),\widehat{e}_{a,k}^{(2)}(X^{(k)}),$ and
$\widehat{p}^{(2)}_{k}(X^{(k)})$. Calculate

$\widehat{\psi}_{\text{\tiny
DR}}^{(1)}(a)=\sum_{k=1}^{K}\widehat{w}_{k}\frac{1}{\sum_{i\in\mathcal{B}_{1}}I(S_{i}=0)}\sum_{i\in\mathcal{B}_{1}}\left(\widehat{o}_{a}^{(2)}(X^{(k)}_{i},A_{i},S_{i})(Y_{i}-\widehat{g}_{a,k}^{(2)}(X^{(k)}_{i}))+I(S_{i}=0)\widehat{g}_{a,k}^{(2)}(X^{(k)}_{i})\right),$

where $\widehat{o}_{a}^{(2)}(X^{(k)},A,S)$ is defined by replacing
$\widehat{e}_{a,k}(X^{(k)})$ and $\widehat{p}_{k}(X^{(k)})$ by
$\widehat{e}_{a,k}^{(2)}(X^{(k)})$ and $\widehat{p}^{(2)}_{k}(X^{(k)})$ in the
definition of $\widehat{o}_{a}(X^{(k)},A,S)$.

  3. 3.

Calculate $\widehat{\psi}_{\text{\tiny DR}}^{(2)}(a)$ by repeating Step 2 with
the roles of $\mathcal{B}_{1}$ and $\mathcal{B}_{2}$ switched.

  4. 4.

The sample split doubly robust estimator is defined as
$\widehat{\psi}_{\text{\tiny
DR}}^{(SS)}(a)=\frac{1}{2}\left(\widehat{\psi}_{\text{\tiny
DR}}^{(1)}(a)+\widehat{\psi}_{\text{\tiny DR}}^{(2)}(a)\right)$.

In Supplementary Web Appendix (A.4) we show that if assumptions B1, B3, and B4
hold (i.e., without requiring the Donsker assumption) the sample split doubly
robust estimator is consistent and has the same asymptotic representation as
in Theorem 2. Hence, the sample split doubly robust estimator is model and
rate doubly robust without requiring any Donsker assumptions.

So far we have discussed the properties the weights
$\\{\widehat{w}_{1},\ldots,\widehat{w}_{k}\\}$ need to satisfy, without
describing a procedure for choosing their values. In traditional individual
participant data meta-analysis, the optimal choice of weights (defined in
terms of minimum asymptotic variance) is to weight each study by the inverse
of the study-specific variance estimator [37]. This suggests using weights
that are proportional to the inverse of the estimated pattern-specific
sampling variance for each term in the sum in equation (9) [38]. The
optimality of inverse variance weights in meta-analysis, however, relies on
independence between the study-specific estimators. In causally interpretable
meta-analysis, trials with each missingness pattern are “standardized” to the
same target population, which induces dependence between the estimators for
different missingness patterns. As a result, inverse variance weighting will
not be in general optimal. In Supplementary Web Appendix A.5 we show that the
weights that minimize the asymptotic variance of the doubly robust estimator
involve solving an estimating equation that includes a term that depends on
the asymptotic variance of the estimator for each missingness pattern and a
term that represents the correlation induced by standardizing the estimators
using the same sample from the target population.

###  4.3 Inference

If for each $k\in\\{1,\ldots,K\\}$, the estimators
$\widehat{g}_{a,k}(X^{(k)})$, $\widehat{p}_{k}(X^{(k)})$, and
$\widehat{e}_{a,k}(X^{(k)})$ converge at a fast enough rate such that
$REM=o_{p}(1)$, then the estimator $\widehat{\psi}_{\text{\tiny DR}}(a)$ is
asymptotically normally distributed. Confidence intervals for
$\widehat{\psi}_{\text{\tiny DR}}(a)$ can be constructed using either sandwich
variance estimators or using the nonparametric bootstrap.

###  4.4 Comparing treatments that are not compared head-to-head in any trial

The methods developed can easily be extended to compare target population
efficacy of treatments that are not compared head-to-head in any trial (i.e.,
to perform “indirect treatment comparisons” as in “network meta-analysis”
[39]). More specifically, we can restrict the estimation of each potential
outcome mean to the set of trials that include the treatment under
consideration. In conventional or pooled meta-analyses, it is difficult to
contrast potential outcome means across studies because differences could
result from differences in treatment efficacy or differences in the underlying
study populations. On the contrary, the methods developed here allow for
comparisons of treatments in the context of the same target population,
provided the identifiability conditions hold.

##  5 Simulations

We conducted simulations to evaluate finite sample performance of the three
estimators proposed in Section 4. We compared the g-formula estimator
$\widehat{\psi}_{{\tiny g}}(a)$, the weighting estimator
$\widehat{\psi}_{\text{\tiny w}}(a)$, and the doubly robust estimator
$\widehat{\psi}_{\text{\tiny DR}}(a)$. For comparison, we also implemented a
naive pooled estimator that simply averages the outcomes for those individuals
assigned to treatment $a$ in the combined data from all trials. That is, the
naive pooled estimator is calculated using the formula
$\widehat{\psi}_{\text{\tiny
pool}}(a)=\frac{\sum_{i=1}^{n}Y_{i}I(A_{i}=a,S_{i}\neq
0)}{\sum_{i=1}^{n}I(A_{i}=a,S_{i}\neq 0)}$. We also implemented “complete-
case” versions of the g-formula, weighting, and the doubly robust estimators
that only used data from the trials that collected information on all
covariates (these estimators ignore all information from trials in which
covariate data are systematically missing).

###  5.1 Data generation

We simulated the covariate vector from a five-dimensional mean-zero
multivariate normal distribution with covariance matrix with element $(i,j)$
equal to $0.6^{|i-j|}$, for $i=1,\ldots,5$, and $j=1,\ldots,5$. Participation
in any trial was simulated from a logistic regression model with

$\Pr[S=0|X]=1-\frac{\exp(1+0.2X_{1}+0.2X_{2}+0.2X_{3}+0.1X_{1}^{2}+0.1X_{2}^{2}+0.1X_{3}^{2})}{1+\exp(1+0.2X_{1}+0.2X_{2}+0.2X_{3}+0.1X_{1}^{2}+0.1X_{2}^{2}+0.1X_{3}^{2})}.$

Observations with $S\neq 0$, were allocated to one of three trials using a
multinomial logistic regression

$S|(X,S\neq
0)\sim\text{Multinomial}\left((p_{1},p_{2},p_{3}),\sum_{i=1}^{n}I(S_{i}\neq
0)\right),$

where

$\displaystyle p_{1}$ $\displaystyle=\Pr[S=1|X,S\neq
0]=\frac{\theta}{1+\theta+\zeta}$ $\displaystyle p_{2}$
$\displaystyle=\Pr[S=2|X,S\neq 0]=\frac{\zeta}{1+\theta+\zeta}$ $\displaystyle
p_{3}$ $\displaystyle=\Pr[S=3|X,S\neq 0]=1-p_{1}-p_{2},$

with

$\displaystyle\theta$
$\displaystyle=\text{exp}\left(\log(1.3)X_{1}+\log(1.3)X_{2}+\log(1.3)X_{3}\right)$
$\displaystyle\zeta$
$\displaystyle=\text{exp}\left(\log(0.8)X_{1}+\log(0.8)X_{2}+\log(0.8)X_{3}\right).$

For all three trials, treatment assignment was simulated from a Bernoulli
distribution with parameter $0.5$, reflecting marginal randomization.

We generated potential outcomes using

$\displaystyle Y^{1}$
$\displaystyle=1+0.2X_{1}+0.2X_{2}+0.1X_{1}^{2}+0.1X_{2}^{2}+\varepsilon^{1}$
$\displaystyle Y^{0}$ $\displaystyle=-0.2X_{1}-0.2X_{4}I(S\neq
2)-0.1X_{1}^{2}-0.1X_{2}^{2}+\varepsilon^{0},$

where $\varepsilon^{a}$ are independent $\mathcal{N}(0,1)$ for
$a\in\\{0,1\\}$. The observed outcome was generated under consistency as
$Y=Y^{1}A+Y^{0}(1-A)$. We considered two systematically missing data
scenarios. In the first missing data scenario, all covariates were fully
observed in trials $S=1$ and $S=2$, but $X_{5}$ was missing in trial $S=3$.
This results in two different missing data patterns ($K=2$). For the second
missing data scenario, all covariates were observed in trial $S=1$, $X_{4}$
was missing in trial $S=2$, and $X_{5}$ was missing in trial $S=3$. This
results in three different missing data patterns ($K=3$). The total sample
size in the collection of trials and the target population was $n=2,000$ with
an average size of the target population sample of $454$ and the average
sample size in trials $1$, $2$, and $3$ of $540$, $540$, and $466$,
respectively. We ran $1000$ simulations for each setting. The potential
outcome means in the target population were
$\operatorname{\textnormal{\mbox{E}}}[Y^{1}|S=0]=1.07$ and
$\operatorname{\textnormal{\mbox{E}}}[Y^{0}|S=0]=-0.094$ (calculated using a
numerical approximation).

###  5.2 Implementation and results

For trial participation, we used (1) a logistic regression model that included
linear and quadratic terms for all covariates that are collected in each trial
and (2) a logistic regression model that only included linear terms. Because
the model for participation in any trial is a logistic regression model that
includes linear and quadratic terms, we expect that, for each missingness-
pattern, a logistic regression model that includes both linear and quadratic
terms can better approximate the correct model whereas the logistic regression
model that includes only linear terms is more severely misspecified. For the
outcome, we used a correctly specified linear regression model that included
linear and quadratic main effects for all covariates that are collected in the
trial and a misspecified model that only included linear terms. As the
treatment assignment was marginally randomized, any model that includes an
intercept is correctly specified. In the simulations we used a main effects
logistic regression estimator to estimate
$\Pr_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]$. Estimating the treatment
assignment can improve precision in randomized controlled trials compared to
using the true value [40].

In the main part of the manuscript, to combine estimates across missing data
patterns, we use weights that are proportional to the total sample size of the
trials in a given missing data pattern $\mathcal{S}^{k}$. Additional
simulations show that selecting the weights proportional to the sample size,
proportional to the inverse of the estimated variance for each term in
equation (9), or using the optimal weights in terms of minimizing the
asymptotic variance resulted in almost identical performance in terms of both
bias and standard deviation (results not shown).

Table 1 shows results from $1,000$ simulations where the estimators are
compared in terms of bias and standard deviation. The naive pooled estimator
was biased for the potential outcome mean in the target population for both
missing data scenarios and both treatments. The estimators
$\widehat{\psi}_{{\tiny g}}(a)$, $\widehat{\psi}_{\text{\tiny w}}(a)$, and
$\widehat{\psi}_{\text{\tiny DR}}(a)$ were approximately unbiased when the
outcome model was correctly specified and the model for trial participation
included both linear and quadratic terms. When the outcome model was
misspecified, the g-formula estimator was biased and when the probability of
trial participation included only first order terms the weighting estimator
was biased. The doubly robust estimator was unbiased if either the outcome
model was correctly specified or the model for trial participation included
both linear and quadratic terms. The g-formula estimator had the smallest
variance, followed by the doubly robust estimator, and the weighting estimator
had the largest variance. Using only data from trials with complete covariate
information resulted in loss in efficiency for all three estimators compared
to using the estimators described in Section 4.1.

$K$ |  Estimator |  Correct Outcome model |  Quadratic Trial Participation Model |  Restricted to Complete Trials | $a=0$ | $a=1$
---|---|---|---|---|---|---
Bias | SD | Bias | SD
2 | Naive |  |  |  | -5.81 | 3.04 | 7.60 | 3.18
GF | x |  |  | -0.44 | 1.83 | 0.094 | 1.97
GF |  |  |  | -1.74 | 1.92 | 1.39 | 1.97
GF | x |  | x | -0.67 | 2.19 | 0.13 | 2.33
W |  | x |  | -0.54 | 1.97 | 0.089 | 1.97
W |  |  |  | -2.64 | 2.01 | 2.32 | 1.96
W |  | x | x | -0.67 | 2.32 | 0.089 | 2.28
DR | x | x |  | -0.44 | 1.88 | 0.085 | 1.97
DR |  | x |  | -0.45 | 1.92 | -0.098 | 1.96
DR | x |  |  | -0.45 | 1.88 | 0.098 | 1.97
DR |  |  |  | -3.00 | 1.96 | 2.60 | 1.93
DR | x | x | x | -0.67 | 2.24 | 0.12 | 2.28
3 | Naive |  |  |  | -5.81 | 3.04 | 7.60 | 3.18
GF | x |  |  | -0.45 | 1.88 | 0.11 | 1.98
GF |  |  |  | -1.74 | 1.92 | 1.39 | 1.97
GF | x |  | x | -0.020 | 3.13 | 0.094 | 3.04
W |  | x |  | -0.58 | 2.01 | 0.094 | 1.96
W |  |  |  | -2.68 | 2.06 | 2.37 | 1.97
W |  | x | x | -0.058 | 3.35 | 0.098 | 3.13
DR | x | x |  | -0.45 | 1.88 | 0.098 | 1.96
DR |  | x |  | -0.44 | 1.92 | 0.11 | 1.98
DR | x |  |  | -0.45 | 1.87 | 0.11 | 1.98
DR |  |  |  | -3.13 | 1.97 | 2.73 | 1.95
DR | x | x | x | -0.020 | 3.18 | 0.089 | 3.09
Table 1: Comparisons of the naive pooled $\widehat{\psi}_{\text{\tiny
pool}}(a)$ (Naive), g-formula $\widehat{\psi}_{{\tiny g}}(a)$ (GF), weighting
$\widehat{\psi}_{\text{\tiny w}}(a)$ (W), and doubly robust
$\widehat{\psi}_{\text{\tiny DR}}(a)$ (DR) estimators in terms of bias and
standard deviation (SD). To facilitate numerical comparisons, the bias and
standard deviation estimators are multiplied by $\sqrt{n}=44.7$. Restricted to
complete trials refers to using only trials with no systematic missing data in
the analysis. Correct outcome model refers to using a correctly specified
outcome model and quadratic trial participation model refers to using a model
for trial participation that includes linear and quadratic terms. $K$ refers
to the the number of missingness-patterns.

##  6 Meta-analysis of lung cancer screening trials

To illustrate the application of the methods, we use data from two lung cancer
screening trials to estimate a potential outcome mean associated with chest
radiography screening in a nationally representative target population of
individuals who would be eligible for lung cancer screening.

###  6.1 Description of trial and target data

The National Lung Screening trial (NLST) enrolled people aged 55 to 74 that
had $\geq$ 30 pack-year history who were current smokers or had quit within
the past 15 years. Participants were randomized to screening with computed
tomography or chest radiography [41, 42]. The trial showed substantial
benefits of lung cancer screening and its results have informed national lung
cancer screening guidelines [43, 44]. The Prostate, Lung, Colorectal, and
Ovarian (PLCO) Cancer Screening Trial enrolled participants aged 55-74 with no
prior history of prostate, lung, colorectal, or ovarian cancer [45, 46]. The
trial randomized participants to either a control arm that received usual care
or a treatment arm where participants received screening for several cancers
including screening for lung cancer using chest radiography. We used the
binary outcome of whether a participant was diagnosed with lung cancer within
six years from study entry.

The target population data came from the 2003-2004 NHANES, a cross-sectional
survey with a multi-stage clustering design evaluating the health and
nutritional status of non-institutionalized US adults and children [47].
NHANES oversampled certain subgroups (including both racial and ethnic
subgroups), resulting in data that require the use of sampling weights to
represent the underlying target population [48]. In the next subsection we
describe how the g-formula, weighting, and doubly robust estimators can be
modified to handle a weighted sampling design in the target population.

The target population we focused on is people recommended for lung cancer
screening in US. As the NLST eligibility criteria were very similar to the
population recommended for screening by the U.S. Preventive Services Task
Force guidelines [43, 44], we restricted the NHANES data to only include
individuals who met the NLST eligibility criteria (and have data on smoking
habits that allow us to verify the eligibility criteria). Table 2 in
Supplementary Web Appendix C shows summary statistics of the covariates used
in the analysis stratified by data-source. The table shows that NHANES
participants are on average less educated, have more comorbidities than
participants in the two trials, and are less likely to live with a smoker
compared to NLST participants. Because both trials evaluated chest radiography
screening, we focused on estimating the potential outcome mean of that
intervention in the NHANES target population (because screening with
computerized tomography and usual care were each evaluated in only one of the
trials, previously described transportability methods [9] would suffice to
transport potential outcome means under these interventions).

Information on whether a participant lived with a smoker and if the
participant had a history of asthma was collected in both the NLST and NHANES
data but not in the PLCO data, resulting in systematically missing data; all
three data sources collected the other covariate information listed in Table 2
in Supplementary Web Appendix C. Furthermore, both NLST and NHANES collected
average number of cigarettes smoked per day as a continuous variable while in
PLCO it was collected as a categorical variable. Thus, the transportability
analysis from NLST to NHANES used average number of cigarettes smoked per day
as a continuous variable while the transportability analysis from PLCO to
NHANES used it as a categorical variable. For simplicity, we restrict the
analysis to observations that have no within trial or target missing covariate
data and trial participants that were not censored in the first six years.
This resulted in $22,841$ observations from NLST, $17,639$ observations from
PLCO, and $219$ observations from NHANES (representing roughly $8.5$ million
participants) being used in our analysis. To implement the estimators, we used
main effects logistic regression for the model for trial participation, the
treatment assignment, and the outcome model, where the model for trial
participation is weighted by the sampling weights as described in the next
section.

###  6.2 Accounting for survey sampling weights and multi-cluster sampling in
the target population

Let $\eta_{i}$ denote the survey sampling weight associated with observation
$i$ in the target population. We set these weights to one for all observations
in the collection of trials. To account for the sampling weights in the data
from the target population we modify the g-formula estimator:

$\widehat{\psi}^{\eta}_{{\tiny
g}}(a)=\sum_{k=1}^{K}\widehat{w}_{k}\frac{1}{\sum_{i=1}^{n}\eta_{i}I(S_{i}=0)}\sum_{i=1}^{n}\eta_{i}I(S_{i}=0)\widehat{g}_{a,k}(X_{i}^{(k)}).$
(11)

Because $\widehat{g}_{a,k}(X^{(k)})$ is an estimator for
$\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},A=a,S\in\mathcal{S}^{k}]$,
it is estimated using only trial data; therefore, the procedure for estimating
$\widehat{g}_{a,k}(X^{(k)})$ does not need to be modified to account for the
sampling weights.

Furthermore, we also modify the weighting estimator to account for the
sampling weights:

$\widehat{\psi}^{\eta}_{\text{\tiny
w}}(a)=\sum_{k=1}^{K}\widehat{w}_{k}\frac{1}{\sum_{i=1}^{n}\eta_{i}I(S_{i}=0)}\sum_{i=1}^{n}\frac{(1-\widetilde{p}_{k}(X^{(k)}_{i}))I(S_{i}\in\mathcal{S}^{k},A_{i}=a)}{\widehat{e}_{a,k}(X_{i}^{(k)})\widetilde{p}_{k}(X^{(k)}_{i})}Y_{i}.$
(12)

Here, the estimator $\widetilde{p}_{k}(X^{(k)})$ for
$\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]$ uses both data from the trials and
target population and needs to account for the weighted sampling design. In
our analysis we used a weighted logistic regression model with weights equal
to one for NLST and PLCO participants and the survey sampling weights for the
NHANES participants. The estimator $\widehat{e}_{a,k}(X^{(k)})$ is an
estimator for $\Pr_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]$ and is only fit using
trial data and does therefore not need to be modified to account for the
sampling weights in the target population data.

Last, we modify the doubly robust estimator as follows:

$\displaystyle\widehat{\psi}^{\eta}_{\text{\tiny
DR}}(a)=\sum_{k=1}^{K}\widehat{w}_{k}$
$\displaystyle\frac{1}{\sum_{i=1}^{n}\eta_{i}I(S_{i}=0)}$
$\displaystyle\times\sum_{i=1}^{n}\Bigg{(}\frac{(1-\widetilde{p}_{k}(X^{(k)}_{i}))I(S_{i}\in\mathcal{S}^{k},A_{i}=a)}{\widehat{e}_{a,k}(X_{i}^{(k)})\widetilde{p}_{k}(X^{(k)}_{i})}(Y_{i}-\widehat{g}_{a,k}(X^{(k)}_{i}))+\eta_{i}I(S_{i}=0)\widehat{g}_{a,k}(X^{(k)}_{i})\Bigg{)}.$

We used the bootstrap for variance estimation and confidence interval (CI)
construction. We accounted for the NHANES complex sampling design by using a
stratified bootstrap [49, 50] where the resampling is done so it is consistent
with the NHANES sampling design (i.e., the resampling is done at the primary
sampling unit and strata level [48]). For the naive pooled estimator, we used
the non-parametric bootstrap stratified by trial.

###  6.3 Results

We estimated the six year risk of being diagnosed with lung cancer to be
$0.073$ (95% bootstrap interval $[0.065,0.082]$) using the g-formula
estimator; $0.072$ (95% bootstrap interval $[0.065,0.079]$) using the
weighting estimator; and $0.072$ (95% bootstrap interval $[0.064,0.080]$)
using the doubly robust estimator. For comparison, the naive pooled estimator
using only the trial data produced an estimated risk of $0.050$ (95% bootstrap
interval $[0.048,0.053]$). Thus, all three estimators for the potential
outcome mean in the target population that we proposed produced very similar
point estimates that were substantially larger compared with the naive pooled
estimator. In fact, the bootstrap interval for the naive pooled estimator did
not overlap with the bootstrap interval for the other three estimators. This
may be due to differences in the distribution of education between the trials
and the target population (see Table 2 in Supplementary Web Appendix C for the
distribution of education levels) and higher education levels being associated
with lower risk of lung cancer diagnosis (see Table 3 in Supplementary Web
Appendix C for prevalence rates by education level).

##  7 Discussion

We provided identifiability results and proposed three estimators for the
potential outcome mean in the target population for causally interpretable
meta-analysis with systematically missing data. We studied the large-sample
properties of the estimators and illustrated that the estimators have good
finite-sample performance in simulation studies. Last, we estimated the risk
of lung cancer diagnosis when chest radiography screening is applied using
data from two large randomized controlled trials and using target population
data from the NHANES study.

Our methods rely on the untestable missing at random assumption (Assumption
A6), it would be of interest to develop sensitivity analysis methods to
evaluate how violations of that assumption may affect the results [51]. We
only considered the setting where the only source of missing data was
systematic missing data and combining the methods we develop here with methods
for handling within-trial or within-target missing data is of interest. Future
work could also consider extensions of the methods to address failure time
outcomes and covariate measurement error.

Throughout, we focused on meta-analyses that aim to estimate causal estimands
such as potential outcome means and average treatment effects in the target
population. In related work, Kundu et al.[52] developed a generalized method
of moments approach for combining information from multiple parametric
regression models to estimate a regression coefficient in the presence of
systematically missing data. Their approach rests on the assumption that the
joint distribution of the outcome and covariates is the same across all the
studies being combined, and requires data from a “reference sample” that can
be used to estimate the joint distribution of the covariates. It would be
interesting to explore whether the approach we used for causal estimands can
be extended to meta-analyses of regression models.

The estimators we proposed combine data from all trials within a given
missingness pattern to estimate the potential outcome mean in the target
population and then form a convex combination of the pattern-specific
estimates to construct a summary estimate. An alternative approach is to
“transport” each trial separately and then use our methods to calculate the
summary estimate combining the trial-specific estimators [38]. If that
approach is taken, the methods we propose here can use robust and efficient
trial-specific estimators [10], combining their estimates using optimal
weights, and accounting for the correlation induced by standardization to the
same target population. In other words, transporting each trial separately can
be viewed as a special case of our methods when each trial collects different
covariates. For instance, this kind of analysis may be particularly attractive
in multi-cohort observational analyses where the contributing cohorts have
different data structures. In addition, working with each trial separately may
be appealing even when some trials have the same missingness pattern when, for
each trial, either the outcome, or both the trial participation and treatment
models can be correctly specified, whereas correctly specifying at least one
group of models for the pooled data from all trials with a given missingness
pattern is more challenging. We note, however, that separately transporting
each trial requires stronger positivity conditions because conditions A1 and
A5 need to hold for each trial rather than the aggregate of all trials with
the same missingness pattern.

## References

  * [1] Michael Borenstein, Larry V Hedges, Julian PT Higgins, and Hannah R Rothstein.  Introduction to meta-analysis.  John Wiley & Sons, 2011.
  * [2] Christopher H Schmid, Theo Stijnen, and Ian White.  Handbook of Meta-Analysis.  CRC Press, 2020.
  * [3] Issa J Dahabreh, Lucia C Petito, Sarah E Robertson, Miguel A Hernán, and Jon A Steingrimsson.  Toward causally interpretable meta-analysis: Transporting inferences from multiple randomized trials to a new target population.  Epidemiology (Cambridge, Mass.), 31(3):334–344, 2020.
  * [4] Issa J Dahabreh and Miguel A Hernán.  Extending inferences from a randomized trial to a target population.  European Journal of Epidemiology, pages 1–4, 2019.
  * [5] Stephen R Cole and Elizabeth A Stuart.  Generalizing evidence from randomized clinical trials to target populations: the ACTG 320 trial.  American Journal of Epidemiology, 172(1):107–115, 2010.
  * [6] Elizabeth Tipton.  Improving generalizations from experiments using propensity score subclassification assumptions, properties, and contexts.  Journal of Educational and Behavioral Statistics, 38(3):239–266, 2012.
  * [7] Colm O’Muircheartaigh and Larry V Hedges.  Generalizing from unrepresentative experiments: a stratified propensity score approach.  Journal of the Royal Statistical Society. Series C (Applied Statistics), 63(2):195–210, 2014.
  * [8] Judea Pearl and Elias Bareinboim.  External validity: from do-calculus to transportability across populations.  Statistical Science, 29(4):579–595, 2014.
  * [9] Issa J Dahabreh, Sarah E Robertson, Jon A Steingrimsson, Elizabeth A Stuart, and Miguel A Hernán.  Extending inferences from a randomized trial to a new target population.  Statistics in Medicine, 39(14):1999–2014, 2020.
  * [10] Issa J Dahabreh, Sarah E Robertson, Lucia C Petito, Miguel A Hernán, and Jon A Steingrimsson.  Efficient and robust methods for causally interpretable meta-analysis: transporting inferences from multiple randomized trials to a target population.  arXiv preprint arXiv:1908.09230, 2019.
  * [11] David H Barker, Issa J Dahabreh, Jon A Steingrimsson, Christopher Houck, Geri Donenberg, Ralph DiClemente, and Larry K Brown.  Causally interpretable meta-analysis: Application in adolescent hiv prevention.  Prevention Science, pages 1–12, 2021.
  * [12] Donald B Rubin.  Inference and missing data.  Biometrika, 63(3):581–592, 1976.
  * [13] Geert Molenberghs, Garrett Fitzmaurice, Michael G Kenward, Anastasios Tsiatis, and Geert Verbeke.  Handbook of missing data methodology.  CRC Press, 2014.
  * [14] 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.
  * [15] Matthieu Resche-Rigon, Ian R White, Jonathan W Bartlett, Sanne AE Peters, Simon G Thompson, and PROG-IMT Study Group.  Multiple imputation for handling systematically missing confounders in meta-analysis of individual participant data.  Statistics in Medicine, 32(28):4890–4905, 2013.
  * [16] Shahab Jolani, Thomas PA Debray, Hendrik Koffijberg, Stef van Buuren, and Karel GM Moons.  Imputation of systematically missing predictors in an individual participant data meta-analysis: a generalized approach using mice.  Statistics in Medicine, 34(11):1841–1863, 2015.
  * [17] Deborah Kunkel and Eloise E Kaizar.  A comparison of existing methods for multiple imputation in individual participant data meta-analysis.  Statistics in Medicine, 36(22):3507–3532, 2017.
  * [18] Shahab Jolani.  Hierarchical imputation of systematically and sporadically missing data: An approximate bayesian approach using chained equations.  Biometrical Journal, 60(2):333–351, 2018.
  * [19] Matthieu Resche-Rigon and Ian R White.  Multiple imputation by chained equations for systematically and sporadically missing multilevel data.  Statistical Methods in Medical Research, 27(6):1634–1649, 2018\.
  * [20] Barry C Arnold and S James Press.  Compatible conditional distributions.  Journal of the American Statistical Association, 84(405):152–156, 1989.
  * [21] Barry C Arnold, Enrique Castillo, and José María Sarabia.  Exact and near compatibility of discrete conditional distributions.  Computational statistics & data analysis, 40(2):231–252, 2002\.
  * [22] Donald B Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  Journal of Educational Psychology, 66(5):688, 1974.
  * [23] James M Robins and Sander Greenland.  Causal inference without counterfactuals: comment.  Journal of the American Statistical Association, 95(450):431–435, 2000.
  * [24] Issa J Dahabreh, James M Robins, Sebastien JP Haneuse, and Miguel A Hernán.  Generalizing causal inferences from randomized trials: counterfactual and graphical identification.  arXiv preprint arXiv:1906.10792, 2019.
  * [25] James M Robins.  A new approach to causal inference in mortality studies with a sustained exposure period – application to control of the healthy worker survivor effect.  Mathematical Modelling, 7(9):1393–1512, 1986.
  * [26] James Robins, Mariela Sued, Quanhong Lei-Gomez, and Andrea Rotnitzky.  Comment: Performance of double-robust estimators when” inverse probability” weights are highly variable.  Statistical Science, 22(4):544–559, 2007.
  * [27] Issa J Dahabreh, Miguel A Hernán, Sarah E Robertson, Ashley Buchanan, and Jon A Steingrimsson.  Generalizing trial findings in nested trial designs with sub-sampling of non-randomized individuals.  arXiv preprint arXiv:1902.06080, 2019.
  * [28] Heejung Bang and James M Robins.  Doubly robust estimation in missing data and causal inference models.  Biometrics, 61(4):962–973, 2005.
  * [29] J Hájek.  Comment on “An essay on the logical foundations of survey sampling by D. Basu”.  In V P Godambe and D A Sprott, editors, Foundations of statistical inference, page 236. Holt, Rinehart, and Winston, New York City, NY, 1971.
  * [30] James M Robins, Mariela Sued, Quanhong Lei-Gomez, and Andrea Rotnitzky.  Comment: Performance of double-robust estimators when “inverse probability” weights are highly variable.  Statistical Science, 22(4):544–559, 2007.
  * [31] Aad W van der Vaart and Jon A Wellner.  Weak Convergence and Empirical Processes.  Springer, 1996.
  * [32] Joel L Horowitz.  Semiparametric and nonparametric methods in econometrics, volume 12.  Springer, 2009.
  * [33] Andrea Rotnitzky, Ezequiel Smucler, and James M Robins.  Characterization of parameters with a mixed bias property.  Biometrika, 108(1):231–238, 2021.
  * [34] 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.
  * [35] 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.
  * [36] 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.
  * [37] D Zeng and DY Lin.  On random-effects meta-analysis.  Biometrika, 102(2):281–294, 2015.
  * [38] Tat-Thang Vo, Raphael Porcher, Anna Chaimani, and Stijn Vansteelandt.  Rethinking meta-analysis: assessing case-mix heterogeneity when combining treatment effects across patient populations.  arXiv preprint arXiv:1908.10613, 2019.
  * [39] Thomas Lumley.  Network meta-analysis for indirect treatment comparisons.  Statistics in Medicine, 21(16):2313–2324, 2002.
  * [40] Jon Arni Steingrimsson, Daniel F Hanley, and Michael Rosenblum.  Improving precision by adjusting for prognostic baseline variables in randomized trials with binary outcomes, without regression model assumptions.  Contemporary Clinical Trials, 54:18–24, 2017.
  * [41] National Lung Screening Trial Research Team.  The national lung screening trial: overview and study design.  Radiology, 258(1):243–253, 2011.
  * [42] National Lung Screening Trial Research Team Writing committee:, Denise R Aberle, Amanda M Adams, Christine D Berg, Jonathan D Clapp, Kathy L Clingan, Ilana F Gareen, David A Lynch, Pamela M Marcus, and Paul F Pinsky.  Baseline characteristics of participants in the randomized national lung screening trial.  Journal of the National Cancer Institute, 102(23):1771–1779, 2010\.
  * [43] Virginia A Moyer.  Screening for lung cancer: Us preventive services task force recommendation statement.  Annals of Internal Medicine, 160(5):330–338, 2014.
  * [44] Alex H Krist, Karina W Davidson, Carol M Mangione, Michael J Barry, Michael Cabana, Aaron B Caughey, et al.  Screening for lung cancer: Us preventive services task force recommendation statement.  JAMA, 325(10):962–970, 2021.
  * [45] Martin M Oken, Willam G Hocking, Paul A Kvale, Gerald L Andriole, Saundra S Buys, Timothy R Church, E David Crawford, Mona N Fouad, Claudine Isaacs, Douglas J Reding, et al.  Screening by chest radiograph and lung cancer mortality: the prostate, lung, colorectal, and ovarian (plco) randomized trial.  JAMA, 306(17):1865–1873, 2011.
  * [46] Martin M Oken, Pamela M Marcus, Ping Hu, Thomas M Beck, William Hocking, Paul A Kvale, Jill Cordes, Thomas L Riley, Stephen D Winslow, Steven Peace, et al.  Baseline chest radiograph for lung cancer detection in the randomized prostate, lung, colorectal and ovarian cancer screening trial.  Journal of the National Cancer Institute, 97(24):1832–1839, 2005\.
  * [47] Clifford L Johnson, Ryne Paulose-Ram, Cynthia L Ogden, Margaret D Carroll, Deanna Kruszan-Moran, Sylvia M Dohrmann, and Lester R Curtin.  National health and nutrition examination survey. analytic guidelines, 1999-2010.  2013\.
  * [48] Lester R Curtin, Leyla K Mohadjer, Sylvia M Dohrmann, Jill M Montaquila, Deanna Kruszan-Moran, Lisa B Mirel, Margaret D Carroll, Rosemarie Hirsch, Susan Schober, and Clifford L Johnson.  The national health and nutrition examination survey: Sample design, 1999-2006.  Vital and health statistics. Series 2, Data evaluation and methods research, (155):1–39, 2012.
  * [49] Jun Shao.  Impact of the bootstrap on sample surveys.  Statistical Science, 18(2):191–198, 2003.
  * [50] Jon NK Rao and CFJ Wu.  Resampling inference with complex survey data.  Journal of the American Statistical Association, 83(401):231–241, 1988.
  * [51] Issa J Dahabreh, James M Robins, Sebastien JP Haneuse, Iman Saeed, Sarah E Robertson, Elisabeth A Stuart, and Miguel A Hernán.  Sensitivity analysis using bias functions for studies extending inferences from a randomized trial to a target population.  arXiv preprint arXiv:1905.10684, 2019.
  * [52] Prosenjit Kundu, Runlong Tang, and Nilanjan Chatterjee.  Generalized meta-analysis for multiple regression models across studies with disparate covariate information.  Biometrika, 106(3):567–585, 2019.
  * [53] Aad W Van der Vaart.  Asymptotic statistics, volume 3.  Cambridge University Press, 2000.

##  Appendix A Proofs

###  A.1 Proof of identifiability results

###### Proof.

Assumptions A1 through A6 imply that for each $k\in\\{1,\ldots,K\\}$ we can
write the potential outcome mean in the target population
$\operatorname{\textnormal{\mbox{E}}}[Y^{a}|S=0]$ as

$\displaystyle\operatorname{\textnormal{\mbox{E}}}[Y^{a}|S=0]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}[\operatorname{\textnormal{\mbox{E}}}[Y^{a}|X,S=0]|S=0]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y^{a}|X,S=0]|S=0]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y^{a}|X,S\in\mathcal{S}^{k}]|S=0]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y^{a}|X,S\in\mathcal{S}^{k},A=a]|S=0]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X,S\in\mathcal{S}^{k},A=a]|S=0]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},S\in\mathcal{S}^{k},A=a]|S=0]$
$\displaystyle=\psi_{k}(a),$

where all expectations above are well defined under positivity conditions A3
and A5.

From the above result, and the constraint $\sum_{k=1}^{K}w_{k}^{*}=1$, it
follows that

$\displaystyle\psi(a)=\sum_{k=1}^{K}w_{k}^{*}\psi_{k}(a)=\sum_{k=1}^{K}w_{k}^{*}\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},A=a,S\in\mathcal{S}^{k}]|S=0]=\operatorname{\textnormal{\mbox{E}}}[Y^{a}|S=0].$

Turning our attention to the weighting re-expression of the identifiability
result, for each $k\in\\{1,\ldots,K\\}$, we have

$\displaystyle\psi_{k}(a)$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},A=a,S\in\mathcal{S}^{k}]|S=0]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}\left[\operatorname{\textnormal{\mbox{E}}}_{k}\left[\frac{YI(A=a,S\in\mathcal{S}^{k})}{\Pr_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]}\Bigg{|}X^{(k)}\right]\Bigg{|}S=0\right]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}\left[\frac{I(S=0)}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\left[\frac{YI(A=a,S\in\mathcal{S}^{k})}{\Pr_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]}\Bigg{|}X^{(k)}\right]\right]$
$\displaystyle=\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\left[\operatorname{\textnormal{\mbox{E}}}_{k}\left[\frac{YI(A=a,S\in\mathcal{S}^{k})(1-\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}])}{\Pr_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]}\Bigg{|}X^{(k)}\right]\right]$
$\displaystyle=\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\left[\frac{YI(A=a,S\in\mathcal{S}^{k})\Pr_{k}[S=0|X^{(k)}]}{\Pr_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]}\right].$

∎

Thus, it follows that

$\psi(a)=\sum_{k=1}^{K}w_{k}^{*}\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\left[\frac{I(A=a,S\in\mathcal{S}^{k})\Pr_{k}[S=0|X^{(k)}]Y}{\Pr_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]}\right].$

###  A.2 First Order Influence Function

We will use pathwise derivatives to calculate the first order influence
function of $\psi(a)$ using the identifiability result in expression (3) (see
e.g., [53] for more details on influcence function calculations). Let $p_{t}$
denote a one dimensional submodel with $t\in[0,1[$ and $t=0$ corresponds to
the true data law. Using that differentiation is a linear operator

$\displaystyle\frac{\partial\psi_{p_{t}}(a)}{\partial
t}\bigg{|}_{t=0}=\sum_{k=1}^{K}w_{k}\frac{\partial}{\partial
t}\operatorname{\textnormal{\mbox{E}}}_{k,p_{t}}[\operatorname{\textnormal{\mbox{E}}}_{k,p_{t}}[Y|X^{(k)},A=a,I(S\in\mathcal{S}^{k})=1]|S=0]\Bigg{|}_{t=0}.$

Result in Appendix D in [10] show that

$\displaystyle\frac{\partial}{\partial
t}\operatorname{\textnormal{\mbox{E}}}_{k,p_{t}}[\operatorname{\textnormal{\mbox{E}}}_{k,p_{t}}[Y|X^{(k)},A=a,I(S\in\mathcal{S}^{k})=1]|S=0]\Bigg{|}_{t=0}$
$\displaystyle\quad=\operatorname{\textnormal{\mbox{E}}}_{k,p_{0}}\Bigg{[}\dfrac{1}{\Pr_{k,p_{0}}[S=0]}\Bigg{\\{}I(S=0)\Big{\\{}\operatorname{\textnormal{\mbox{E}}}_{k,p_{0}}[Y|X^{(k)},I(S\in\mathcal{S}^{k})=1,A=a]-\psi_{p_{0}}(a)\Big{\\}}$
$\displaystyle\quad\quad\quad\quad+\dfrac{I(S\in\mathcal{S}^{k},A=a)\Pr_{k,p_{0}}[S=0|X^{(k)}]}{\Pr_{k,p_{0}}[S\in\mathcal{S}^{k}|X^{(k)}]\Pr_{k,p_{0}}[A=a|X^{k},I(S\in\mathcal{S}^{k})=1]}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\times\Big{\\{}Y-\operatorname{\textnormal{\mbox{E}}}_{k,p_{0}}[Y|X^{(k)},I(S\in\mathcal{S}^{k})=1,A=a]\Big{\\}}\Bigg{\\}}u(O)\Bigg{]},$

where $u(O)$ is the score of the observable data. It follows that for any set
of weights $\\{w_{1},\ldots,w_{K}\\}$ the influence function of $\psi(a)$
under the non-parametric model is

$\displaystyle IF=$
$\displaystyle\sum_{k=1}^{K}w_{k}\dfrac{1}{\Pr_{k,p_{0}}[S=0]}\Bigg{\\{}I(S=0)\Big{\\{}\operatorname{\textnormal{\mbox{E}}}_{k,p_{0}}[Y|X^{(k)},I(S\in\mathcal{S}^{k})=1,A=a]-\psi_{p_{0}}(a)\Big{\\}}$
$\displaystyle\quad+\dfrac{I(S\in\mathcal{S}^{k},A=a)\Pr_{k,p_{0}}[S=0|X^{(k)}]}{\Pr_{k,p_{0}}[S\in\mathcal{S}^{k}|X^{(k)}]\Pr_{k,p_{0}}[A=a|X^{k},I(S\in\mathcal{S}^{k})=1]}$
$\displaystyle\quad\quad\quad\times\Big{\\{}Y-\operatorname{\textnormal{\mbox{E}}}_{k,p_{0}}[Y|X^{(k)},I(S\in\mathcal{S}^{k})=1,A=a]\Big{\\}}\Bigg{\\}}.$
(A.1)

###  A.3 Proof of Theorem 2

As $n\rightarrow\infty$ we have

$\displaystyle\widehat{\psi}_{\text{\tiny DR}}(a)\overset{P}{\longrightarrow}$
$\displaystyle\sum_{k=1}^{K}w^{*}_{k}\Bigg{(}\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\Bigg{[}\frac{(1-p^{*}_{k}(X^{(k)}))I(S\in\mathcal{S}^{k},A=a)}{e_{a,k}^{*}(X^{(k)})p^{*}_{k}(X^{(k)})}\big{(}Y-g_{a,k}^{*}(X^{(k)})\big{)}$
$\displaystyle\quad\quad\quad\quad\quad+I(S=0)g_{a,k}^{*}(X^{(k)})\Bigg{]}\Bigg{)}.$

We will now show that the right hand side of the above equation is equal to
$\psi(a)$ for the two cases listed in Assumption B1.

Case 1: Let $k\in\\{1,\ldots,K\\}$ be given and assume
$\widehat{g}_{a,k}(X^{(k)})\overset{P}{\longrightarrow}\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},S\in\mathcal{S}^{k},A=a]$.
We don’t assume that either of the estimators $\widehat{p}_{k}(X^{(k)})$ or
$\widehat{e}_{a,k}(X^{(k)})$ are consistent. Under this assumptions we have

$\displaystyle\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\left[I(S=0)g^{*}_{a}(X^{(k)})\right]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},S\in\mathcal{S}^{k},A=a]|S=0]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}[Y^{a}|S=0],$

and

$\displaystyle\operatorname{\textnormal{\mbox{E}}}_{k}\Bigg{[}\frac{(1-p^{*}_{k}(X^{(k)}))I(S\in\mathcal{S}^{k},A=a)}{e_{a,k}^{*}(X^{(k)})p^{*}_{k}(X^{(k)})}\big{(}Y-g_{a,k}^{*}(X^{(k)})\big{)}\Bigg{]}$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}\Bigg{[}\operatorname{\textnormal{\mbox{E}}}_{k}\Bigg{[}\frac{(1-p^{*}_{k}(X^{(k)}))I(S\in\mathcal{S}^{k},A=a)}{e_{a,k}^{*}(X^{(k)})p^{*}_{k}(X^{(k)})}\big{(}Y-\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},S\in\mathcal{S}^{k},A=a]\big{)}\Bigg{|}X^{(k)}\Bigg{]}\Bigg{]}$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}\Bigg{[}\frac{(1-p^{*}_{k}(X^{(k)}))\Pr_{k}[S\in\mathcal{S}^{k},A=a|X^{(k)}]}{e_{a,k}^{*}(X^{(k)})p^{*}_{k}(X^{(k)})}\big{(}Y-\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},S\in\mathcal{S}^{k},A=a]\big{)}\Bigg{|}X^{(k)},S\in\mathcal{S}^{k},A=a\Bigg{]}$
$\displaystyle=0$

Case 2: Let $k\in\\{1,\ldots,K\\}$ be given, and assume
$\widehat{p}_{k}(X^{(k)})\overset{P}{\longrightarrow}\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]$
and
$\widehat{e}_{a,k}(X^{(k)})\overset{P}{\longrightarrow}\Pr_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]$.
We don’t assume that $\widehat{g}_{a,k}(X^{(k)})$ is a consistent estimator.
Under this assumption

$\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\left[I(S=0)g^{*}_{a}(X^{(k)})\right]=\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\left[{\Pr}_{k}[S=0|X^{(k)}]g^{*}_{a}(X^{(k)})\right]$

and

$\displaystyle\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\Bigg{[}$
$\displaystyle\frac{(1-p^{*}_{k}(X^{(k)}))I(S\in\mathcal{S}^{k},A=a)}{e_{a,k}^{*}(X^{(k)})p^{*}_{k}(X^{(k)})}g_{a,k}^{*}(X^{(k)})\Bigg{]}$
$\displaystyle=\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\Bigg{[}\frac{(1-p^{*}_{k}(X^{(k)}))g_{a,k}^{*}(X^{(k)})}{e_{a,k}^{*}(X^{(k)})p^{*}_{k}(X^{(k)})}E_{k}[I(S\in\mathcal{S}^{k},A=a)|X^{(k)}]\Bigg{]}$
$\displaystyle=\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\left[{\Pr}_{k}[S=0|X^{(k)}]g^{*}_{a}(X^{(k)})\right].$

We have

$\displaystyle\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\Bigg{[}$
$\displaystyle\frac{(1-p^{*}_{k}(X^{(k)}))I(S\in\mathcal{S}^{k},A=a)}{e_{a,k}^{*}(X^{(k)})p^{*}_{k}(X^{(k)})}Y\Bigg{]}$
$\displaystyle=\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\big{[}\Pr_{k}[S=0|X^{(k)}]\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},A=a,S\in\mathcal{S}^{k}]\big{]}$
$\displaystyle=\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\big{[}I(S=0)\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},A=a,S\in\mathcal{S}^{k}]\big{]}$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},A=a,S\in\mathcal{S}^{k}]|S=0]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}[Y^{a}|S=0]$
$\displaystyle=\psi(a).$

Combing the above, we have shown that for both cases and each
$k\in\\{1,\ldots,K\\}$,

$\frac{1}{\Pr_{k}[S=0]}\operatorname{\textnormal{\mbox{E}}}_{k}\Bigg{[}\frac{(1-p^{*}_{k}(X^{(k)}))I(S\in\mathcal{S}^{k},A=a)}{e_{a,k}^{*}(X^{(k)})p^{*}_{k}(X^{(k)})}\big{(}Y-g_{a,k}^{*}(X^{(k)})\big{)}+I(S=0)g_{a,k}^{*}(X^{(k)})\Bigg{]}=\psi(a).$

As the weights sum to one we get

$\widehat{\psi}_{\text{\tiny
DR}}(a)\overset{P}{\longrightarrow}\sum_{k=1}^{K}w^{*}_{k}\psi(a)=\psi(a),$

completing the consistency proof.

Now consider the asymptotic representation given by equation (10). Rewrite

$\displaystyle\sqrt{n}\big{(}\widehat{\psi}_{\text{\tiny
DR}}(a)-\psi(a)\big{)}$
$\displaystyle=\sum_{k=1}^{K}w_{k}^{*}\mathbb{G}_{k}\bigg{(}L(\widehat{g}_{a,k}(X^{(k)}),\widehat{e}_{a,k}(X^{(k)}),\widehat{\gamma}_{k},\widehat{p}_{k}(X^{(k)}))\bigg{)}$
$\displaystyle-\sum_{k=1}^{K}w_{k}^{*}\mathbb{G}_{k}\bigg{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}_{k}(X^{(k)}))\bigg{)}$
$\displaystyle+\sum_{k=1}^{K}w_{k}^{*}\mathbb{G}_{k}\bigg{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}_{k}(X^{(k)}))\bigg{)}$
$\displaystyle+\sqrt{n}\left(\sum_{k=1}^{K}w_{k}^{*}\operatorname{\textnormal{\mbox{E}}}\bigg{[}L(\widehat{g}_{a,k}(X^{(k)}),\widehat{e}_{a,k}(X^{(k)}),\widehat{\gamma}_{k},\widehat{p}_{k}(X^{(k)}))\bigg{]}-\sum_{k=1}^{K}w_{k}^{*}\psi(a)\right)$
$\displaystyle+\sqrt{n}\sum_{k=1}^{K}(\widehat{w}_{k}-w^{*}_{k})\mathbb{P}_{k}\big{(}L(\widehat{g}_{a,k}(X^{(k)}),\widehat{e}_{a,k}(X^{(k)}),\widehat{\gamma}_{k},\widehat{p}_{k}(X^{(k)}))\big{)}$

By the Donsker condition

$\displaystyle\mathbb{G}_{k}\bigg{(}$ $\displaystyle
L(\widehat{g}_{a,k}(X^{(k)}),\widehat{e}_{a,k}(X^{(k)}),\widehat{\gamma}_{k},\widehat{p}_{k}(X^{(k)}))\bigg{)}$
$\displaystyle-\mathbb{G}_{k}\bigg{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}_{k}(X^{(k)}))\bigg{)}=o_{p}(1).$

Similar calculations to Appendix E in [10] show that for each
$k\in\\{1,\ldots,K\\}$,

$\displaystyle E\big{[}w^{*}_{k}$ $\displaystyle
L(\widehat{g}_{a,k}(X^{(k)}),\widehat{e}_{a,k}(X^{(k)}),\widehat{\gamma}_{k},\widehat{p}_{k}(X^{(k)}))\big{]}-w_{k}^{*}\psi(a)$
$\displaystyle\leq
O_{P}\Bigg{(}\Big{(}\Big{|}\Big{|}\big{(}\widehat{p}_{k}(X^{(k)})-{\Pr}_{k}[S\in\mathcal{S}^{k}|X^{(k)}]\Big{|}\Big{|}_{2}+\Big{|}\Big{|}\widehat{e}_{a,k}(X^{(k)})-{\Pr}_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]\Big{|}\Big{|}_{2}\Big{)}$
$\displaystyle\times\Big{|}\Big{|}\widehat{g}_{a,k}(X^{(k)})-\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},S\in\mathcal{S}^{k},A=a]\Big{|}\Big{|}_{2}\Bigg{)},$
(A.2)

which completes the proof of the asymptotic representation.

###  A.4 Asymptotic representation of sample split doubly robust estimator

As the consistency proof of Theorem 2 does not require the Donsker condition,
we only need to show that the asymptotic representation provided in Theorem 2
holds for the sample split doubly robust estimator without requiring the
Donsker assumption.

Recall that sample split doubly robust estimator is given by

$\widehat{\psi}_{\text{\tiny
DR}}^{(SS)}(a)=\frac{1}{2}\widehat{\psi}_{\text{\tiny
DR}}^{(1)}(a)+\frac{1}{2}\widehat{\psi}_{\text{\tiny DR}}^{(2)}(a),$

where

$\widehat{\psi}_{\text{\tiny
DR}}^{(1)}(a)=\sum_{k=1}^{K}\widehat{w}_{k}\frac{1}{\sum_{i\in\mathcal{B}_{1}}I(S_{i}=0)}\sum_{i\in\mathcal{B}_{1}}\left(\widehat{o}_{a}^{(2)}(X^{(k)}_{i},A_{i},S_{i})(Y_{i}-\widehat{g}_{a,k}^{(2)}(X^{(k)}_{i}))+I(S_{i}=0)\widehat{g}_{a,k}^{(2)}(X^{(k)}_{i})\right)$

and

$\widehat{\psi}_{\text{\tiny
DR}}^{(2)}(a)=\sum_{k=1}^{K}\widehat{w}_{k}\frac{1}{\sum_{i\in\mathcal{B}_{2}}I(S_{i}=0)}\sum_{i\in\mathcal{B}_{2}}\left(\widehat{o}_{a}^{(1)}(X^{(k)}_{i},A_{i},S_{i})(Y_{i}-\widehat{g}_{a,k}^{(1)}(X^{(k)}_{i}))+I(S_{i}=0)\widehat{g}_{a,k}^{(1)}(X^{(k)}_{i})\right).$

For simplicity assume both set of indexes $\mathcal{B}_{1}$ and
$\mathcal{B}_{2}$ are of equal size $m=\frac{n}{2}$. Define
$\mathbb{P}_{k}^{(1)}(W)$ and $\mathbb{P}_{k}^{(2)}(W)$ as the analogs of
$\mathbb{P}_{k}(W)$ where the average is taken over observations in
$\mathcal{B}_{1}$ and $\mathcal{B}_{2}$, respectively. Define
$\mathbb{G}_{k}^{(1)}(W)=\sqrt{m}(\mathbb{P}_{k}^{(1)}(W)-E[W])$ and
$\mathbb{G}_{k}^{(2)}(W)=\sqrt{m}(\mathbb{P}_{k}^{(2)}(W)-E[W])$. For $j=1,2$
and $k\in\\{1,\ldots,K\\}$, let
$\gamma_{k}^{(j)}=\Pr_{k}[S=0|\mathcal{B}_{j}]^{-1}$,
$\widehat{\gamma}_{k}^{(1)}=\frac{\sum_{i\in\mathcal{B}_{1}}(I(S_{i}=0)+I(S_{i}\in\mathcal{S}^{k}))}{\sum_{i\in\mathcal{B}_{1}}I(S_{i}=0)}$
and
$\widehat{\gamma}_{k}^{(2)}=\frac{\sum_{i\in\mathcal{B}_{2}}(I(S_{i}=0)+I(S_{i}\in\mathcal{S}^{k}))}{\sum_{i\in\mathcal{B}_{2}}I(S_{i}=0)}$.
Rewrite

$\displaystyle\sqrt{n}\left(\widehat{\psi}_{\text{\tiny
DR}}^{(SS)}(a)-\psi(a)\right)$
$\displaystyle=\frac{\sqrt{n}}{2}\left(\widehat{\psi}_{\text{\tiny
DR}}^{(1)}(a)-\psi(a)\right)+\frac{\sqrt{n}}{2}\left(\widehat{\psi}_{\text{\tiny
DR}}^{(2)}(a)-\psi(a)\right)$
$\displaystyle=\frac{\sqrt{n}}{2}\sum_{k=1}^{K}(\widehat{w}_{k}-w^{*}_{k})\mathbb{P}^{(1)}_{n,k}\big{(}L(\widehat{g}_{a,k}^{(2)}(X^{(k)}),\widehat{e}_{a,k}^{(2)}(X^{(k)}),\widehat{\gamma}_{k}^{(2)},\widehat{p}^{(2)}_{k}(X^{(k)}))\big{)}$
$\displaystyle+\frac{\sqrt{n}}{2}\sum_{k=1}^{K}(\widehat{w}_{k}-w^{*}_{k})\mathbb{P}^{(2)}_{n,k}\big{(}L(\widehat{g}_{a,k}^{(1)}(X^{(k)}),\widehat{e}_{a,k}^{(1)}(X^{(k)}),\widehat{\gamma}_{k}^{(1)},\widehat{p}^{(1)}_{k}(X^{(k)}))\big{)}$
$\displaystyle+\frac{\sqrt{n}}{2}\left(\sum_{k=1}^{K}w_{k}^{*}\mathbb{P}^{(1)}_{n,k}\big{(}L(\widehat{g}_{a,k}^{(2)}(X^{(k)}),\widehat{e}_{a,k}^{(2)}(X^{(k)}),\widehat{\gamma}_{k}^{(2)},\widehat{p}^{(2)}_{k}(X^{(k)}))\big{)}-\psi(a)\right)$
$\displaystyle+\frac{\sqrt{n}}{2}\left(\sum_{k=1}^{K}w_{k}^{*}\mathbb{P}^{(2)}_{n,k}\big{(}L(\widehat{g}_{a,k}^{(1)}(X^{(k)}),\widehat{e}_{a,k}^{(1)}(X^{(k)}),\widehat{\gamma}_{k}^{(1)},\widehat{p}^{(1)}_{k}(X^{(k)}))\big{)}-\psi(a)\right)$
$\displaystyle=\frac{1}{\sqrt{2}}\Bigg{(}\sum_{k=1}^{K}w_{k}^{*}\bigg{(}\mathbb{G}^{(1)}_{n,k}\big{(}L(\widehat{g}_{a,k}^{(2)}(X^{(k)}),\widehat{e}_{a,k}^{(2)}(X^{(k)}),\widehat{\gamma}_{k}^{(2)},\widehat{p}^{(2)}_{k}(X^{(k)}))\big{)}$
$\displaystyle-\mathbb{G}^{(1)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}\bigg{)}$
$\displaystyle+\sum_{k=1}^{K}w_{k}^{*}\mathbb{G}^{(1)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}$
$\displaystyle+\sqrt{m}\left(\operatorname{\textnormal{\mbox{E}}}\left[\sum_{k=1}^{K}w_{k}^{*}L(\widehat{g}_{a,k}^{(2)}(X^{(k)}),\widehat{e}_{a,k}^{(2)}(X^{(k)}),\widehat{\gamma}_{k}^{(2)},\widehat{p}^{(2)}_{k}(X^{(k)}))\right]-\psi(a)\right)\Bigg{)}$
$\displaystyle+\frac{1}{\sqrt{2}}\Bigg{(}\sum_{k=1}^{K}w_{k}^{*}\bigg{(}\mathbb{G}^{(2)}_{n,k}\big{(}L(\widehat{g}_{a,k}^{(1)}(X^{(k)}),\widehat{e}_{a,k}^{(1)}(X^{(k)}),\widehat{\gamma}_{k}^{(1)},\widehat{p}^{(1)}_{k}(X^{(k)}))\big{)}$
$\displaystyle-\mathbb{G}^{(2)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}\bigg{)}$
$\displaystyle+\sum_{k=1}^{K}w_{k}^{*}\mathbb{G}^{(2)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}$
$\displaystyle+\sqrt{m}\left(\operatorname{\textnormal{\mbox{E}}}\left[\sum_{k=1}^{K}w_{k}^{*}L(\widehat{g}_{a,k}^{(1)}(X^{(k)}),\widehat{e}_{a,k}^{(1)}(X^{(k)}),\widehat{\gamma}_{k}^{(1)},\widehat{p}^{(1)}_{k}(X^{(k)}))\right]-\psi(a)\right)\Bigg{)}$
$\displaystyle+O_{P}\left(\sqrt{n}\sum_{k=1}^{K}(\widehat{w}_{k}-w_{k}^{*})\right).$

By rearranging terms we get

$\displaystyle\sqrt{n}\left(\widehat{\psi}_{\text{\tiny
DR}}^{(SS)}(a)-\psi(a)\right)$
$\displaystyle=\frac{1}{\sqrt{2}}\Big{(}\sum_{k=1}^{K}w_{k}^{*}\mathbb{G}^{(1)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}$
$\displaystyle\quad\quad+\sum_{k=1}^{K}w_{k}^{*}\mathbb{G}^{(2)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}\Big{)}$
$\displaystyle+\frac{1}{\sqrt{2}}\Bigg{(}\sum_{k=1}^{K}w_{k}^{*}\bigg{(}\mathbb{G}^{(1)}_{n,k}\big{(}L(\widehat{g}_{a,k}^{(2)}(X^{(k)}),\widehat{e}_{a,k}^{(2)}(X^{(k)}),\widehat{\gamma}_{k}^{(2)},\widehat{p}^{(2)}_{k}(X^{(k)}))\big{)}$
$\displaystyle-\mathbb{G}^{(1)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}$
$\displaystyle+\sum_{k=1}^{K}w_{k}^{*}\bigg{(}\mathbb{G}^{(2)}_{n,k}\big{(}L(\widehat{g}_{a,k}^{(1)}(X^{(k)}),\widehat{e}_{a,k}^{(1)}(X^{(k)}),\widehat{\gamma}_{k}^{(1)},\widehat{p}^{(1)}_{k}(X^{(k)}))\big{)}$
$\displaystyle-\mathbb{G}^{(2)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}\bigg{)}\Bigg{)}$
$\displaystyle+\frac{1}{\sqrt{2}}\sqrt{m}\left(\operatorname{\textnormal{\mbox{E}}}\left[\sum_{k=1}^{K}w_{k}^{*}L(\widehat{g}_{a,k}^{(2)}(X^{(k)}),\widehat{e}_{a,k}^{(2)}(X^{(k)}),\widehat{\gamma}_{k}^{(1)},\widehat{p}^{(2)}_{k}(X^{(k)}))\right]-\psi(a)\right)$
$\displaystyle+\frac{1}{\sqrt{2}}\sqrt{m}\left(\operatorname{\textnormal{\mbox{E}}}\left[\sum_{k=1}^{K}w_{k}^{*}L(\widehat{g}_{a,k}^{(1)}(X^{(k)}),\widehat{e}_{a,k}^{(1)}(X^{(k)}),\widehat{\gamma}_{k}^{(1)},\widehat{p}^{(1)}_{k}(X^{(k)}))\right]-\psi(a)\right)$
$\displaystyle+O_{P}\left(\sqrt{n}\sum_{k=1}^{K}(\widehat{w}_{k}-w_{k}^{*})\right).$

We have

$\displaystyle\mathbb{P}^{(1)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}+\mathbb{P}^{(2)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}$
$\displaystyle=\frac{2}{n}\sum_{i\in\mathcal{B}_{1}}L(g_{a,k}^{*}(X_{i}^{(k)}),e_{a,k}^{*}(X_{i}^{(k)}),\gamma_{k},p^{*}(X_{i}^{(k)}))+\frac{2}{n}\sum_{i\in\mathcal{B}_{2}}L(g_{a,k}^{*}(X_{i}^{(k)}),e_{a,k}^{*}(X_{i}^{(k)}),\gamma_{k},p^{*}(X_{i}^{(k)}))$
$\displaystyle=2\mathbb{P}_{k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}.$

Using the result above

$\displaystyle\frac{1}{\sqrt{2}}\Big{(}\mathbb{P}^{(1)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}-\operatorname{\textnormal{\mbox{E}}}\left[L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\right]$
$\displaystyle+\mathbb{P}^{(2)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}-\operatorname{\textnormal{\mbox{E}}}\left[L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\right]\Big{)}$
$\displaystyle=\sqrt{2}\mathbb{P}_{k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}-\sqrt{2}\operatorname{\textnormal{\mbox{E}}}\left[L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\right].$

It follows that

$\displaystyle\frac{1}{\sqrt{2}}\left(\sum_{k=1}^{K}w_{k}^{*}\mathbb{G}^{(1)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}+\sum_{k=1}^{K}w_{k}^{*}\mathbb{G}^{(2)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}\right)$
$\displaystyle=\frac{\sqrt{m}}{\sqrt{2}}\Bigg{(}\sum_{k=1}^{K}w_{k}^{*}\mathbb{P}^{(1)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}-\sum_{k=1}^{K}w_{k}^{*}\operatorname{\textnormal{\mbox{E}}}\left[L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\right]\Bigg{)}$
$\displaystyle+\frac{\sqrt{m}}{\sqrt{2}}\Bigg{(}\sum_{k=1}^{K}w_{k}^{*}\mathbb{P}^{(2)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}-\sum_{k=1}^{K}w_{k}^{*}\operatorname{\textnormal{\mbox{E}}}\left[L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\right]\Bigg{)}$
$\displaystyle=\sqrt{n}\left(\sum_{k=1}^{K}w_{k}^{*}\mathbb{P}_{k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}-\sum_{k=1}^{K}w_{k}^{*}\operatorname{\textnormal{\mbox{E}}}\left[L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\right]\right)$
$\displaystyle=\sum_{k=1}^{K}w_{k}^{*}\mathbb{G}_{k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}.$

For each $k\in\\{1,\ldots,K\\}$ the inequality in (A.2) and that $n=2m$ gives

$\displaystyle\frac{1}{\sqrt{2}}\sqrt{m}\left(\operatorname{\textnormal{\mbox{E}}}\left[\mathbb{P}^{(1)}_{n,k}\big{(}L(\widehat{g}_{a,k}^{(2)}(X^{(k)}),\widehat{e}_{a,k}^{(2)}(X^{(k)}),\widehat{\gamma}_{k}^{(1)},\widehat{p}^{(2)}_{k}(X^{(k)}))\big{)}\right]-\psi(a)\right)$
$\displaystyle+\frac{1}{\sqrt{2}}\sqrt{m}\left(\operatorname{\textnormal{\mbox{E}}}\left[\mathbb{P}^{(2)}_{n,k}\big{(}L(\widehat{g}_{a,k}^{(1)}(X^{(k)}),\widehat{e}_{a,k}^{(1)}(X^{(k)}),\widehat{\gamma}_{k}^{(1)},\widehat{p}^{(1)}_{k}(X^{(k)}))\big{)}\right]-\psi(a)\right)$
$\displaystyle\leq\sqrt{n}O_{P}\Bigg{(}\Big{(}\Big{|}\Big{|}\big{(}\widehat{p}_{k}(X^{(k)})-\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]\Big{|}\Big{|}_{2}+\Big{|}\Big{|}\widehat{e}_{a,k}(X^{(k)})-\Pr_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]\Big{|}\Big{|}_{2}\Big{)}$
$\displaystyle\times\Big{|}\Big{|}\widehat{g}_{a,k}(X^{(k)})-\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},S\in\mathcal{S}^{k},A=a]\Big{|}\Big{|}_{2}\Bigg{)}.$

For each $k\in\\{1,\ldots,K\\}$ define

$\displaystyle C_{1,k}=$
$\displaystyle\mathbb{G}^{(1)}_{n,k}\big{(}L(\widehat{g}_{a,k}^{(2)}(X^{(k)}),\widehat{e}_{a,k}^{(2)}(X^{(k)}),\widehat{\gamma}_{k}^{(2)},\widehat{p}^{(2)}_{k}(X^{(k)}))\big{)}-\mathbb{G}^{(1)}_{n,k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}$
$\displaystyle=\sqrt{m}\frac{1}{m}\sum_{i\in\mathcal{B}_{1}}\Bigg{(}L(\widehat{g}_{a,k}^{(2)}(X_{i}^{(k)}),\widehat{e}_{a,k}^{(2)}(X_{i}^{(k)}),\widehat{\gamma}_{k}^{(2)},\widehat{p}^{(2)}_{k}(X_{i}^{(k)}))-L(g_{a,k}^{*}(X_{i}^{(k)}),e_{a,k}^{*}(X_{i}^{(k)}),\gamma_{k},p^{*}_{k}(X_{i}^{(k)}))$
$\displaystyle-\left(\operatorname{\textnormal{\mbox{E}}}\left[L(\widehat{g}_{a,k}^{(2)}(X_{i}^{(k)}),\widehat{e}_{a,k}^{(2)}(X_{i}^{(k)}),\widehat{\gamma}_{k}^{(2)},\widehat{p}^{(2)}_{k}(X_{i}^{(k)}))\right]-\operatorname{\textnormal{\mbox{E}}}\left[L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}_{k}(X^{(k)}))\right]\right)\Bigg{)}$

Define $\mathcal{O}_{1}=\\{O_{i},i\in\mathcal{B}_{1}\\}$ and
$\mathcal{O}_{2}=\\{O_{i},i\in\mathcal{B}_{2}\\}$. As $\mathcal{O}_{1}$ and
$\mathcal{O}_{2}$ are independent, for any $i\in\mathcal{B}_{1}$

$\displaystyle\operatorname{\textnormal{\mbox{E}}}[L(\widehat{g}_{a,k}^{(2)}(X_{i}^{(k)}),\widehat{e}_{a,k}^{(2)}(X_{i}^{(k)}),\widehat{\gamma}_{k}^{(2)},\widehat{p}^{(2)}_{k}(X_{i}^{(k)}))-L(g_{a,k}^{*}(X_{i}^{(k)}),e_{a,k}^{*}(X_{i}^{(k)}),\gamma_{k},p^{*}_{k}(X_{i}^{(k)}))|\mathcal{O}_{2}]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}[L(\widehat{g}_{a,k}^{(2)}(X_{i}^{(k)}),\widehat{e}_{a,k}^{(2)}(X_{i}^{(k)}),\widehat{\gamma}_{k}^{(2)},\widehat{p}^{(2)}_{k}(X_{i}^{(k)}))-L(g_{a,k}^{*}(X_{i}^{(k)}),e_{a,k}^{*}(X_{i}^{(k)}),\gamma_{k},p^{*}_{k}(X_{i}^{(k)}))]$

Hence $\operatorname{\textnormal{\mbox{E}}}[C_{1,k}|\mathcal{O}_{2}]=0.$
Again, using that $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ are independent we
have

$\displaystyle\text{Var}[C_{1,k}|\mathcal{O}_{2}]$
$\displaystyle=\text{Var}\left[L(\widehat{g}_{a,k}^{(2)}(X_{i}^{(k)}),\widehat{e}_{a,k}^{(2)}(X_{i}^{(k)}),\widehat{\gamma}_{k}^{(2)},\widehat{p}^{(2)}_{k}(X_{i}^{(k)}))-L(g_{a,k}^{*}(X_{i}^{(k)}),e_{a,k}^{*}(X_{i}^{(k)}),\gamma_{k},p^{*}_{k}(X_{i}^{(k)}))\bigg{|}\mathcal{O}_{2}\right],$

where $i$ is used to index a random observation with index in
$\mathcal{B}_{1}$. Using this

$\displaystyle\text{Var}[C_{1,k}|\mathcal{O}_{2}]$
$\displaystyle\leq\Pr\left[\left(L(\widehat{g}_{a,k}^{(2)}(X_{i}^{(k)}),\widehat{e}_{a,k}^{(2)}(X_{i}^{(k)}),\widehat{\gamma}_{k}^{(2)},\widehat{p}^{(2)}_{k}(X_{i}^{(k)}))-L(g_{a,k}^{*}(X_{i}^{(k)}),e_{a,k}^{*}(X_{i}^{(k)}),\gamma_{k},p^{*}_{k}(X_{i}^{(k)}))\right)^{2}\Bigg{|}\mathcal{O}_{2}\right].$

By assumptions B3,
$\text{Var}[C_{1,k}|\mathcal{O}_{2}]=\Pr[C_{1,k}^{2}|\mathcal{O}_{2}]\rightarrow
0$ and by Chebyshev’s inequality for any $\delta>0$

$\displaystyle\Pr[|C_{1,k}|>\delta|\mathcal{O}_{2}]\leq\frac{\Pr[C_{1,k}^{2}|\mathcal{O}_{2}]}{\delta^{2}}\rightarrow
0$

when $m\rightarrow\infty$. It follows that
$\sum_{k=1}^{K}w_{k}^{*}C_{1,k}=o_{p}(1)$ and the same arguments with the
roles of $\mathcal{B}_{1}$ and $\mathcal{B}_{2}$ switched shows that
$\sum_{k=1}^{K}w_{k}^{*}C_{2,k}=o_{p}(1)$. Combining all the above gives

$\displaystyle\sqrt{n}\left(\widehat{\psi}_{\text{\tiny
DR}}^{(SS)}(a)-\psi(a)\right)=\sum_{k=1}^{K}w_{k}^{*}\mathbb{G}_{k}\big{(}L(g_{a,k}^{*}(X^{(k)}),e_{a,k}^{*}(X^{(k)}),\gamma_{k},p^{*}(X^{(k)}))\big{)}+REM+o_{p}(1),$

where

$\displaystyle REM\leq$
$\displaystyle\sqrt{n}O_{P}\Bigg{(}\sum_{k=1}^{K}\Big{|}\Big{|}w_{k}^{*}-w^{*}_{k}\Big{|}\Big{|}_{2}+\sum_{k=1}^{K}\Big{|}\Big{|}\widehat{g}_{a,k}(X^{(k)})-\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},S\in\mathcal{S}^{k},A=a]\Big{|}\Big{|}_{2}$
$\displaystyle\times\Big{(}\Big{|}\Big{|}\widehat{p}_{k}(X^{(k)})-{\Pr}_{k}[S\in\mathcal{S}^{k}|X^{(k)}]\Big{|}\Big{|}_{2}+\Big{|}\Big{|}\widehat{e}_{a,k}(X^{(k)})-{\Pr}_{k}[A=a|X^{(k)},S\in\mathcal{S}^{k}]\Big{|}\Big{|}_{2}\Big{)}\Bigg{)}.$

###  A.5 Derivation of optimal weights

All expectations in this section are w.r.t. the biased sampling model
described in Section 3. The expected square of the influence function in
expression (A.1) can be written as

$\displaystyle E[IF^{2};w_{1},\ldots,w_{k}]$
$\displaystyle=\sum_{k=1}^{K}w_{k}^{2}\dfrac{1}{\Pr_{k}[S=0]^{2}}\operatorname{\textnormal{\mbox{E}}}_{k}\Bigg{[}\Bigg{(}I(S=0)\Big{\\{}\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},I(S\in\mathcal{S}^{k})=1,A=a]-\psi(a)\Big{\\}}$
$\displaystyle\quad+\dfrac{I(S\in\mathcal{S}^{k},A=a)\Pr_{k}[S=0|X^{(k)}]}{\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]\Pr_{k}[A=a|X^{k},I(S\in\mathcal{S}^{k})=1]}$
$\displaystyle\quad\quad\quad\times\Big{\\{}Y-\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},I(S\in\mathcal{S}^{k})=1,A=a]\Big{\\}}\Bigg{)}^{2}\Bigg{]}$
$\displaystyle+2\sum_{k=1}^{K}\sum_{j=1}^{k-1}w_{k}w_{j}\dfrac{1}{\Pr_{k}[S=0]\Pr_{j}[S=0]}\operatorname{\textnormal{\mbox{E}}}\big{[}(\operatorname{\textnormal{\mbox{E}}}[Y|X^{(k)},I(S\in\mathcal{S}^{k})=1,A=a]-\psi(a))$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\times
I(S=0)(\operatorname{\textnormal{\mbox{E}}}[Y|X^{(j)},I(S\in\mathcal{S}^{j})=1,A=a]-\psi(a))\big{]}.$
(A.3)

We use the notation
$\operatorname{\textnormal{\mbox{E}}}[IF^{2};w_{1},\ldots,w_{k}]$ to emphasize
the dependence of the influence function on the weights $w_{1},\ldots,w_{k}$.
Define

$\displaystyle V_{k}:=$
$\displaystyle\dfrac{1}{\Pr_{k}[S=0]^{2}}\operatorname{\textnormal{\mbox{E}}}_{k}\Bigg{[}\Bigg{(}I(S=0)\Big{\\{}\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},I(S\in\mathcal{S}^{k})=1,A=a]-\psi(a)\Big{\\}}$
$\displaystyle\quad+\dfrac{I(S\in\mathcal{S}^{k},A=a)\Pr_{k}[S=0|X^{(k)}]}{\Pr_{k}[S\in\mathcal{S}^{k}|X^{(k)}]\Pr_{k}[A=a|X^{k},I(S\in\mathcal{S}^{k})=1]}$
$\displaystyle\quad\quad\quad\times\Big{\\{}Y-\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},I(S\in\mathcal{S}^{k})=1,A=a]\Big{\\}}\Bigg{\\}}\Bigg{)}^{2}\Bigg{]}$

and

$\displaystyle C_{kj}$
$\displaystyle=\dfrac{1}{\Pr_{k}[S=0]\Pr_{j}[S=0]}\operatorname{\textnormal{\mbox{E}}}\big{[}(\operatorname{\textnormal{\mbox{E}}}[Y|X^{(k)},I(S\in\mathcal{S}^{k})=1,A=a]-\psi(a))$
$\displaystyle\quad\quad\quad\quad\times
I(S=0)(\operatorname{\textnormal{\mbox{E}}}[Y|X^{(j)},I(S\in\mathcal{S}^{j})=1,A=a]-\psi(a)\big{]}.$

Using this notation rewrite

$E[IF^{2};w_{1},\ldots,w_{k}]=\sum_{k=1}^{K}w_{k}^{2}V_{k}+2\sum_{k=1}^{K}\sum_{j=1}^{k-1}w_{k}w_{j}C_{kj}.$

To find the optimal weights, defined in terms of minimizing the asymptotic
variance, we need to minimize
$\operatorname{\textnormal{\mbox{E}}}[IF^{2};w_{1},\ldots,w_{k}]$ as a
function of $(w_{1},\ldots,w_{k})$ subject to the constraint
$\sum_{k=1}^{K}w_{k}=1$. Using a Lagrange multiplier $\lambda$ we rewrite the
optimization problem and minimize

$T(w_{1},\ldots,w_{k},\lambda)=\operatorname{\textnormal{\mbox{E}}}[IF^{2};w_{1},\ldots,w_{k}]-\lambda\left(\sum_{k=1}^{K}w_{k}-1\right).$
(A.4)

Or equivalently as

$T(w_{1},\ldots,w_{k},\lambda)=\sum_{k=1}^{K}w_{k}^{2}V_{k}+2\sum_{k=1}^{K}\sum_{j<k}w_{k}w_{k}C_{kj}-\lambda\left(\sum_{k=1}^{K}w_{k}-1\right).$

Both $V_{k}$ and $C_{kj}$ depend on unknown quantities that can be estimated
using plug-in estimators. Denote the estimators by $\widehat{V}_{k}$ and
$\widehat{C}_{kj}$. The empirical version of the Lagrance multiplier equation
is given by

$\widehat{T}(w_{1},\ldots,w_{k},\lambda)=\sum_{k=1}^{K}w_{k}^{2}\widehat{V}_{k}+2\sum_{k=1}^{K}\sum_{j<k}w_{k}w_{j}\widehat{C}_{kj}-\lambda\left(\sum_{k=1}^{K}w_{k}-1\right).$

Differentiating $\widehat{T}(w_{1},\ldots,w_{k},\lambda)$ w.r.t. $w_{j}$, for
$j\in\\{1,\ldots,K\\}$, gives

$\frac{\partial\widehat{T}(w_{1},\ldots,w_{k},\lambda)}{\partial
w_{j}}=2w_{j}\widehat{V}_{j}+2\sum_{k\neq j}w_{k}\widehat{C}_{kj}-\lambda.$
(A.5)

Differentiating $\widehat{T}(w_{1},\ldots,w_{k},\lambda)$ w.r.t. $\lambda$
gives

$-\left(\sum_{k=1}^{K}w_{k}-1\right).$ (A.6)

The optimal weights are obtained by finding the zero crossing of expressions
(A.5) and (A.6).

##  Appendix B Sampling model

Let $D$ be the indicator of whether an individual contributes data to the
analyses. Let $s^{k}_{1},\ldots,s^{k}_{M_{k}}$ be the trials in
$\mathcal{S}^{k}$. By the assumptions made

$\operatorname{\textnormal{\mbox{E}}}[Y|X^{(k)},S=s^{k}_{1},A=a]=\ldots=\operatorname{\textnormal{\mbox{E}}}[Y|X^{(k)},S=s^{k}_{M_{k}},A=a]=\operatorname{\textnormal{\mbox{E}}}[Y|X^{(k)},S\in\mathcal{S}^{k},A=a].$

By the random sampling from the population underlying each trial and from the
target population we have

$\displaystyle\psi_{k}(a)$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},A=a,S\in\mathcal{S}^{k}]|S=0]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},A=a,S\in\mathcal{S}^{k},D=1]|S=0,D=1].$

Working with the above expression gives

$\displaystyle\psi_{k}(a)$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}[\operatorname{\textnormal{\mbox{E}}}_{k}[Y|X^{(k)},A=a,S\in\mathcal{S}^{k},D=1]|S=0,D=1]$
$\displaystyle=\operatorname{\textnormal{\mbox{E}}}_{k}\left[\operatorname{\textnormal{\mbox{E}}}_{k}\left[\frac{YI(A=a,S\in\mathcal{S}^{k})}{\Pr[A=a|X^{(k)},S\in\mathcal{S}^{k},D=1]\Pr[S\in\mathcal{S}^{k}|X^{(k)},D=1]}\Bigg{|}X^{(k)},D=1\right]\Bigg{|}S=0,D=1\right]$
$\displaystyle=\frac{1}{\Pr[S=0|D=1]}\operatorname{\textnormal{\mbox{E}}}_{k}\left[\operatorname{\textnormal{\mbox{E}}}_{k}\left[\frac{Y(1-\Pr[S\in\mathcal{S}^{k}|X^{(k)},D=1])I(A=a,S\in\mathcal{S}^{k})}{\Pr[A=a|X^{(k)},S\in\mathcal{S}^{k},D=1]\Pr[S\in\mathcal{S}^{k}|X^{(k)},D=1]}\Bigg{|}X^{(k)},D=1\right]\Bigg{|}D=1\right]$
$\displaystyle=\frac{1}{\Pr[S=0|D=1]}\operatorname{\textnormal{\mbox{E}}}_{k}\left[\frac{Y(1-\Pr[S\in\mathcal{S}^{k}|X^{(k)},D=1])I(A=a,S\in\mathcal{S}^{k})}{\Pr[A=a|X^{(k)},S\in\mathcal{S}^{k},D=1]\Pr[S\in\mathcal{S}^{k}|X^{(k)},D=1]}\Bigg{|}D=1\right]$

All quantities in the above expressions condition on the available data
($D=1$). Thus, the potential outcome mean $\psi(a)$ are identifiable under the
stratified sampling model.

##  Appendix C Additional data analysis results

Table 2 shows the distribution of the covariates in the National Lung
Screening Trial (NLST), Prostate, Lung, Colorectal, and Ovarian trial (PLCO),
and NHANES datasets (weighted for the NHANES data). Table 3 shows the
prevalence rate of lung cancer diagnosis within six years from study
enrollment for the NLST, PLCO trials by education level.

Table 2: Summary of participant characteristics in the National Lung Screening Trial (NLST), the Prostate, Lung, Colorectal, and Ovarian trial (PLCO), and the NHANES data (using the sampling probabilities). Variable | NLST | PLCO | NHANES
---|---|---|---
Age | 61.3 (5.0) | 62.4 (5.2) | 63.1 (5.5)
BMI | 27.9 (5.1) | 27.6 (4.8) | 28.6 (5.6)
Race (White) | 90.1% | 91.0% | 84.5%
Education (Some college education) | 55.5% | 53.2% | 45.7%
Education (High school graduate) | 38.8% | 37.5% | 29.1%
Smoke years | 39.6 (7.3) | 36.1 (9.2) | 42.5 (7.4)
Gender (Male) | 58.7% | 63.8% | 63.0%
Marital status (Married) | 68.4% | 72.8% | 64.7%
Pack year | 55.5 (23.8) | 57.0 (26.2) | 60.6 (28.9)
History of diabetes (Yes) | 9.3% | 8.5% | 20.3%
History of emphysema (Yes) | 7.3% | 6.5% | 8.7%
| History of heart disease
---
or heart attack (Yes)
12.2% | 13.2% | 28.0%
History of hypertension (Yes) | 34.5% | 35.3% | 45.8%
History of asthma (Yes) | 6.2% | NA | 10.7%
Lived with a smoker (Yes) | 87.5% | NA | 49.3%
Cigarettes per day continuous | 28.4 (11.4) |  | 28.7 (12.7)
Cigarettes per day categorical (1-10) |  | 0.03% |
Cigarettes per day categorical (11-20) |  | 36.5% |
Cigarettes per day categorical (21-30) |  | 30.7% |
Cigarettes per day categorical (31-40) |  | 19.4% |
Cigarettes per day categorical (41-60) |  | 10.9% |
Cigarettes per day categorical (61-80) |  | 2.1% |
Cigarettes per day categorical ($>80$ ) |  | 0.4% |
For both trials, the characteristics are restricted to the chest radiography arms. Continuous covariates are summarized as mean (standard deviation) and categorical covariates are summarized as percentage in each category. BMI is body max index; Smoke years is the total number of years the participant smoked cigarettes; Smoke age is the age at smoking onset; Pack years is calculated as (Total number of years Smoked $\times$ Cigarettes Per Day$/20$). NA indicates that the covariate was not collected. Cigarettes per day were collected as a continuous variable in NHANES and NLST and categorical in PLCO. The summaries for the NHANES data are weighted using the NHANES sampling weights. Table 3: Prevalence rate of lung cancer diagnosis within six years from study enrollment for the National Lung Screening Trial (NLST) and the Prostate, Lung, Colorectal, and Ovarian trial (PLCO) trials by education level. | NLST | PLCO
---|---|---
Did not finish high school | 0.049 | 0.010
High school degree | 0.036 | 0.079
Some college education | 0.027 | 0.065