# Inference procedures in sequential trial emulation with survival outcomes:
comparison of confidence intervals based on the sandwich variance estimator
and two forms of bootstrap

Juliette M. Limozin
MRC Biostatistics Unit
University of Cambridge
<EMAIL_ADDRESS>
& Shaun R. Seaman
MRC Biostatistics Unit
University of Cambridge
<EMAIL_ADDRESS>
& Li Su
MRC Biostatistics Unit
University of Cambridge
<EMAIL_ADDRESS>

###### Abstract

The target trial emulation (TTE) framework is a useful approach for providing
causal inference of treatment effects from observational data. In TTE, inverse
probability weighting is often used to adjust for time-varying confounding
and/or dependent censoring. Then structural models for potential outcomes can
be fitted to the weighted observed data to obtain per-protocol treatment
effect estimates. For inference, the sandwich variance estimator with fixed
weights is a popular but often conservative estimator for the variance of the
structural model parameter estimators. Nonparametric bootstrap accounts for
variation in weight estimation, but is computationally expensive. As an
alternative, the linearised estimating function (LEF) bootstrap is
computationally more efficient, but hasn’t been adapted to the TTE settings.
In this paper, we extend the LEF bootstrap to construct confidence intervals
(CIs) for the marginal risk difference of a survival outcome in sequentially
emulated trials. The performance of CIs based on the LEF bootstrap is compared
with those CIs based on nonparametric bootstrap, two types of CIs based on
Jackknife resampling and CIs constructed using the sandwich variance estimator
in various TTE settings through simulations. Our results suggest that on data
with small/moderate sample sizes, the LEF bootstrapping-based CIs had better
coverage whereas sandwich variance estimator-based CIs’ coverage deteriorated
at later follow-ups and had high failure rates. LEF bootstrap-based CIs were
less sensitive to treatment group imbalance compared to nonparametric
bootstrap-based ones, and the former method was on average 3 times faster to
compute than the latter. For data with large sample sizes, all methods amount
to CIs with similar coverage, with the sandwich variance estimator approach
being the fastest to compute. Overall our findings suggest that LEF
bootstrapping can be advantageous for more accurate inference in the context
of survival analysis using sequential trial emulation.

_Keywords_ Causal Inference, Target Trial Emulation, Survival Analysis,
Longitudinal Data, Inverse Probability Weighting, Marginal Structural Models,
Sequential Trial Emulation, Sandwich Variance Estimation, Bootstrap

##  1 Introduction

###  1.1 Target trial emulation with survival outcomes

Target Trial Emulation (TTE) hernan_using_2016 has become a popular approach
for causal inference using observational longitudinal data
hansford_reporting_2023. The goal of TTE is to estimate and make inferences
about causal effects of treatments that are comparable to those that would be
obtained from a target Randomised Controlled Trial (RCT) hernan_using_2016;
matthews_comparing_2021. TTE can be helpful when it is not possible to conduct
this RCT because of time, budget and ethical restrictions. Hernán and Robins
(2016) have proposed a formal framework for TTE hernan_using_2016, which
highlights the need to specify the target trial’s protocol, i.e., the protocol
of the RCT that would have ideally been conducted, in order to guide the
design and analysis of the emulated trial using data extracted from
observational databases, e.g., disease registries or electronic health
records. The protocol should include the eligibility criteria, treatment
strategies, assignment procedures, outcome(s) of interest, follow-up periods,
causal contrast of interest and an analysis plan; see some step-by-step guides
to TTE in hernan_target_2022, matthews_target_2022, maringe_reflection_2020.

There are various sources of bias that must be addressed in TTE. Firstly,
unlike in an RCT, non-random assignment of treatment at baseline must be
accounted for when estimating the causal effect (intention-to-treat or per-
protocol) of a treatment in an emulated trial. Secondly, similar to an RCT, it
is necessary to account for censoring caused by loss to follow-up in an
emulated trial. Thirdly, when the per-protocol effect is the causal effect of
interest, it is also necessary to handle non-adherence to assigned treatments
induced by time-varying confounders. To address these issues in TTE with
survival outcomes, a useful approach is to fit a marginal structural Cox model
using inverse probability weighting, after first artificially censoring the
patients’ follow-up at the time of treatment non-adherence
daniel_methods_2013; hernan_marginal_2000; hernan_marginal_2001;
murray_causal_2021; seaman_review_2013; clare_causal_2019;
hernan_observational_2008; robins_marginal_2000; cain_inverse_2009;
robins_correcting_2000. Baseline confounders are included in this Cox model as
covariates to adjust for the non-random treatment assignment at baseline using
regression. The inverse probability weights are the product of two sets of
time-varying weights: one to address selection bias from censoring due to loss
to follow-up, and one to address selection bias from the artificial censoring
due to treatment non-adherence. Both sets of weights can be estimated using
regression models that include time-varying variables as covariates.
Counterfactual hazard ratios can be estimated in this marginal structural Cox
model. A modification of this approach is to discretise the survival time and
replace the Cox model with a sequence of logistic regression models as the
Marginal Structural Model (MSM) singer_its_1993; hernan_marginal_2000;
hernan_marginal_2001. Provided that the probability of failure between the
discrete times is small, this sequence of logistic models well approximates
the Cox model dagostino_relation_1990.

The counterfactual hazard ratio has been criticised as lacking causal
interpretations and it has been proposed that other estimands be used instead,
e.g., the marginal risk difference (MRD)buchanan_worth_2014;
hernan_hazards_2010; keogh_causal_2023. The MRD over time can be estimated by
first using the hazard ratio estimates from the fitted MSM to predict the
survival probabilities of the patients in the emulated trial under two
scenarios: when all are treated and when none are treated. For each of these
two scenarios, the predicted survival probabilities are averaged over all
enrolled patients. The estimate of the MRD is then calculated as the
difference between these two averages murray_causal_2021.

###  1.2 Constructing confidence intervals in sequential trial emulation

A potential problem when emulating a trial for the purpose of estimating a
treatment effect is that the number of patients eligible for inclusion in the
trial at any given time may be small. This can be addressed by using the
sequential trials approach(hernan_observational_2008; gran_sequential_2010),
which takes advantage of the fact that patients may meet the eligibility
criteria for the target trial multiple times during their follow-up in the
observational database. In this approach, a sequence of target trials is
emulated, each one starting at a different time. The data from these multiple
trials are pooled and then analysed to produce an overall estimate of the
treatment effect (keogh_causal_2023). This approach was first proposed by
Hernán et al. (2008) hernan_observational_2008 and Gran et al. (2010)
gran_sequential_2010 as a simple way to improve the efficiency of treatment
effect estimation relative to emulating only a single trial. There have been
several applications of the sequential trial approach; see
danaei_observational_2013 and an extended list of works reviewed by Keogh et
al. (2023) keogh_causal_2023.

Despite the increasing popularity of the TTE framework, there is a lack of
research on the performance of different methods for constructing confidence
intervals (CIs) for treatment effects when using the sequential trials
approach. To account for the correlations induced by the same patients being
eligible for multiple trials, a variance estimate for the hazard ratio can be
obtained by using the sandwich variance estimator lin_robust_1989 or using
bootstrap. In causal survival analysis with propensity score weighting to
adjust for baseline confounding of point treatments (austin_variance_2016;
shu_variance_2021), the robust sandwich variance estimator proposed by Lin and
Wei lin_robust_1989 is frequently used. However, this variance estimator does
not account for the uncertainty due to the estimation of the weights, and can
consequently overestimate the true variance (austin_variance_2016;
shu_variance_2021). More complex sandwich variance estimators that account for
this uncertainty are available. Specifically, Shu et al. (2021)
shu_variance_2021 proposed a variance estimator for the hazard ratio of a
point treatment in a marginal structural Cox model when inverse probability
weighting (IPW) is used to deal with baseline confounding. However, it is
unclear how their variance estimator can be extended to the sequential trial
emulation setting with IPW for both treatment non-adherence and dependent
censoring. Enders et al. (2018) enders_robust_2018 described a variance
estimator for the hazard ratio parameters of a marginal structural Cox model
when IPW is used to deal with censoring due to loss to follow-up and treatment
switching due to time-varying confounding. They found no substantial
differences in the performance of the robust sandwich variance estimator and
the consistent variance estimator, and the consistent variance estimator
performed comparatively poorly under scenarios with small samples and many
confounders. No off-the-shelf software has implemented this consistent
variance estimator by Enders et al. (2018).

As an alternative approach, bootstrap resampling has been recommended due to
its ability to account for uncertainty in estimating the weights.
Specifically, in simple settings with IPW estimators of point treatment
effects on a survival outcome, Austin (2016) found that bootstrap-based CIs
performed better than the CIs based on the sandwich variance estimator
austin_variance_2016 with moderate sample sizes (1000). For continuous and
binary outcomes with IPW estimators of point treatment effects, Austin (2022)
confirmed that when sample sizes were small (250 or 500) to moderate (1000),
bootstrap resulted in more accurate estimates of standard errors than
asymptotic variance estimators. However, they also found that bootstrap-based
CIs did not achieve nominal coverage when the sample sizes were small to
moderate and the prevalence of treatment was either very low or very high. Mao
et al. (2018) mao_propensity_2018 observed similar results when constructing
CIs of hazard ratios for binary point treatment effect estimation in survival
analysis using IPW: with small (500) and moderate (1000) sample sizes and
strong associations between confounding variables and treatment assignment,
both the sandwich variance estimation and bootstrap variance estimation
methods resulted in CI coverage under the nominal level. For longitudinal
settings with binary time-varying treatments, Seaman and Keogh (2023)
seaman_simulating_2023 found that with moderate (1000) sample size, bootstrap
and the sandwich variance estimator both led to slightly under-coverage of CIs
for hazard ratios in a marginal structural Cox model with IPW but with no
notable difference between the two methods. With small (250, 500) sample
sizes, coverage of the CIs deteriorated, but with bootstrap CI coverage being
closer to the nominal level than the sandwich variance estimator-based CI.

The works mentioned above all focused on the variance estimation and CIs of
counterfactual hazard ratios, which are often chosen as the estimand in the
TTE literature with survival outcomes hernan_observational_2008;
danaei_observational_2013. While these works were not researched specifically
in the context of emulating a target trial for treatment effect estimation
with survival outcomes, they could be easily adapted to such a setting.
Unfortunately however, regardless of the causal estimands for survival
outcomes, less attention appears to be paid to the development and evaluation
of CI construction methods for sequential trial emulation. Specifically, there
is a lack of research on comparing the performance of nonparametric bootstrap
and the sandwich variance estimator for constructing CIs in various settings
of sequential trial emulation in terms of the number of patients who ever
received treatments (i.e., treatment prevalence), strength of association
between time-varying confounders and treatment/outcome (i.e., strength of
time-varying confounding) and outcome event rates. It is also desirable to
develop computationally more efficient CI methods than nonparametric
bootstrap, which is computationally demanding.

###  1.3 Our contribution

To fill this gap, in this paper, we aim to investigate the methods for
constructing CIs of the MRD in sequential trial emulation with survival
outcomes by carrying out an extensive simulation study to evaluate the
performance of four CIs methods. The first method uses the robust sandwich
variance estimator that ignores the uncertainty in estimation of the inverse
probability weights. The second method is nonparametric bootstrap. This has
the drawback of being computationally expensive. For this reason, we also
investigate a computationally less intensive Linearised Estimating Function
(LEF) bootstrap procedure. Hu and Kalbfleisch (2000) hu_estimating_2000 first
developed the estimating function bootstrap approach for computationally
faster implementation of bootstrap. Further, in settings with cross-
sectional/longitudinal survey data with design weights, Rao and Tausi
rao_estimating_2004 and Binder et al. binder_design-based_2004 proposed the
LEF bootstrap to improve computational efficiency and to avoid ill-conditioned
matrices when fitting logistic models to bootstrap samples. We develop two
forms of the LEF bootstrap for constructing CIs of MRD in sequential trial
emulation. The fourth method is Jackknife resampling
tibshirani_introduction_1994; lipsitz_using_1990; friedl_jackknife_2001, which
had been previously used for constructing CIs of the counterfactual hazard
ratio of treatment effect in the sequential trial emulation setting by Gran et
al. (2010) gran_sequential_2010. Jackknife resampling could be advantageous
for small sample sizes as in such cases it is computationally faster to
implement than bootstrap resampling. In our simulation study, we consider
scenarios with varying strengths of time-varying confounding, varying
treatment prevalence, varying outcome event rates and also varying sample
sizes. Our results provide some guidance to practitioners on which CI
construction methods should be favoured in certain settings of sequential
trial emulation with survival outcomes.

The paper is organised as follows. In Section 2 we introduce the HIV
Epidemiology Research Study (HERS) data as a motivating example and describe a
sequential trial emulation protocol based on the HERS data. Section 3
introduces the notation, causal estimand and causal assumptions and also
describes the estimation procedure for MRD. In Section 4, we describe our
proposed LEF bootstrap methods and other CIs construction methods that we
compare in the setting of this paper. In Section 5, we conduct a simulation
study to compare these methods. In Section 6, the methods are applied to the
HERS data. We conclude in Section 7 with a discussion.

##  2 HIV Epidemiology Research Study: a motivating example

This research is motivated by the data from the HIV Epidemiology Research
Study (HERS) ko_estimating_2003, a natural history study that included 1310
women with, or at high risk of, HIV infection at four sites (Baltimore,
Detroit, New York, Providence) enrolled between 1993-1995 and followed up to
2000. The HERS had 12 scheduled visits, one approximately every six months,
where a variety of clinical, behavioural, and sociological outcomes as well as
self-reported treatment information were recorded.

Following Ko et al. (2003) ko_estimating_2003 and Yiu and Su (2022)
yiu_joint_2022, we aim to evaluate the effect of self-reported Highly Active
AntiRetroviral Treatment (HAART) among HIV-infected patients in the HERS
cohort on their risk of death from any cause. HAART is a combination of three
or more antiretroviral therapies that treat HIV infection, and it helps reduce
the number of copies the virus makes in the body to slow down the development
of AIDS eggleton_highly_2022. Other clinical and demographic variables related
to treatment assignment and disease progression were recorded, including CD4
cell count, HIV viral load, self-reported HIV symptoms, race, and the site in
which a patient was enrolled. Following Yiu and Su (2022) yiu_joint_2022, we
treat visit 8 in 1996 as the baseline of the observational cohort, as HAART
was more widely used and recorded in the HERS by then. In total, there were
584 women assessed at visit 8. Their exact time of death during the follow-up
of the HERS was recorded. There were also dropouts due to patients’ loss to
follow-up but they are treated as discrete since exact times are unknown. Yiu
and Su (2022) conducted their analysis with standard MSMs with IPW by defining
the time-varying treatment as ordinal with 3 levels: “no treatment”,
“antiretroviral therapy other than HAART", and “HAART”yiu_joint_2022. In our
analysis we applied the TTE framework with a focus on a binary treatment
strategy of “HAART" versus "no HAART" .

Table 1 outlines a hypothetical RCT that could have been conducted to estimate
the per-protocol effect of HAART (vs. other or no treatment) on all-cause
mortality. This target trial could be emulated using the HERS data and the
sequential trials method outlined in Section 1 and described in more detail in
the next section.

Table 1: Protocol of a target trial to estimate the per-protocol effect of HAART (vs. other or no treatment) on all-cause mortality, which could be emulated using the HERS data Eligibility criteria |  HIV-infected women with no prior history of HAART.
---|---
Treatment strategies |  Patients will be assigned to either the treatment or control group at trial baseline and expected to remain in their assigned group. Those in the treatment group will receive HAART and those in the control group will not receive HAART treatment but may take other antiviral treatments.
Assignment procedures |  Patients will be randomly assigned to either treatment strategy at baseline and will be aware of the strategy to which they have been assigned.
Follow-up period |  Patients are followed from treatment randomisation at baseline visit, for a course of 5 visits spaced 6 months apart, or until event or loss to follow-up, whichever occurs first.
Outcome |  Death from any cause
Causal contrast of interest |  Per-protocol effect of sustained HAART treatment on all-cause mortality
Statistical Methods |  Patients’ follow-up will be artificially censored at the time they deviate from their assigned treatment strategy at baseline. IPW will be used to address the selection bias of this artificial censoring due to treatment non-adherence and confounding by time-varying covariates as well as the selection bias of dependent censoring due to patients’ dropouts.

As mentioned in Section 1, a practical problem for TTE using observational
data is that if we only emulate one trial with a single enrolment period, the
number of patients who initiated the treatment at the baseline and the number
of outcome events among them could be small, which prevents any meaningful
analysis. In the HERS data, only 76 patients are in the HAART treatment arm at
cohort baseline defined at visit 8, versus 415 patients in the control arm,
after applying the eligibility criteria (see Table 1). As mentioned
previously, patients could meet trial eligibility criteria multiple times
during their follow-up in the observational database. In our example, an
additional 62 women initiated the HAART treatment at the first follow-up visit
after cohort baseline. Therefore, it is desirable to emulate a sequence of
target trials and combine their analyses for efficient estimation of treatment
effects.

We could emulate 5 trials labelled from 0 to 4 with sequential enrollment
periods so that the trials start at visits 8, 9, 10, 11 and 12 respectively.
The trial protocol, and more specifically eligibility criteria remain the same
across all 5 trials, which in our example means patients must have no prior
HAART history before the start of each trial to be eligible for enrollment in
either treatment or control arm. The study horizons differ: trial 0, starting
at visit 8, has 5 follow-ups from visit 8 to 12; trial 1 starting at visit 9,
has 4 follow-ups from visit 9 to 12, and so on. This approach means we can use
data from patients who started receiving the HAART later in the observational
cohort: 76 women were eligible for trial 0 in the treated arm, 62 women in
trial 1, 50 women in trial 2, 26 women in trial 3 and 20 women in trial 4. We
have therefore increased the total number of HAART initiators from 76 patients
to 234 patients by using the sequential trial emulation approach. Patients can
be eligible for multiple trials while they still meet the eligibility
criteria. For example, a patient who had not been receiving HAART at visits 8
and 9 but started to receive HAART from visit 10 will be eligible as a member
of the control arm in trials 0 and 1, and as a member of the treatment arm in
trial 2. As data from the same patients can be included in multiple trials,
the correlation between observed data within the same patients will need to be
taken into account in the analysis. Figure 7 in Section A of the Supplementary
Material provides a schematic illustration of the sequential trial emulation
approach.

##  3 Estimation of the per-protocol effect in sequentially emulated trials

###  3.1 Setting and notation

Consider an observational study in which $n$ patients were followed up from
time $t_{0}$ until the earliest of the event of interest, loss to follow-up,
and the end of the study. For each patient, time-independent variables are
measured at time $t_{0}$, and time-varying variables are measured at regular
times $t_{0}<t_{1}<\ldots<t_{n_{v}-1}$ during follow-up, where $t_{n_{v}}$
denotes the time of the end of the study ($t_{{n_{v}}-1}<t_{n_{v}}$). Data on
each patient are assumed to be independent and identically distributed. Data
from this study will be used to create data for a sequence of ${n_{v}}$
trials. The $m$th sequential trial (i.e., trial $m$, $m=0,\ldots,{n_{v}}-1$)
begins at time $t_{m}$, includes patients who are eligible for enrolment at
this time, and ends at time $t_{n_{v}}$. Hence, within this trial, the time-
varying variables are measured at times $t_{m},\ldots,t_{{n_{v}}-1}$. We shall
refer to these ${n_{v}}-m$ measurement times as the trial visits for trial
$m$. For example, trial visit 0 takes place at time $t_{m}$ and trial visit
${n_{v}}-m-1$ takes place at time $t_{{n_{v}}-1}$. For each of the $n$
patients in the observational study, we define the following variables.

  * •

$E_{m}$ is an indicator of whether the patient is eligible ($E_{m}=1$) or not
($E_{m}=0$) for trial $m$.

  * •

$Y_{m,k}=1$ if the patient experiences the event of interest in time interval
$[t_{m},t_{m+k+1})$, and $Y_{m,k}=0$ otherwise.

  * •

$V$ denotes the patient’s vector of time-independent covariates measured at
time $t_{0}$.

  * •

$A_{m,k}$ and $L_{m,k}$ denote, respectively, the patient’s treatment and
time-varying variables measured at time $t_{m+k}$.

  * •

$C_{m,k}$ is an indicator that the patient is censored due to loss to follow-
up in the interval $[t_{m+k},t_{m+k+1})$. So, if $C_{m,k}=1$ then
$Y_{m,k+1},Y_{m,k+2},...,Y_{m+1,k},Y_{m+1,k+1},...$ are not observed.

$Y_{m,k}$, $A_{m,k}$, $L_{m,k}$ and $C_{m,k}$ will serve as, respectively, the
outcome indicator, the binary treatment indicator, the time-varying variables,
and the censoring indicator, at trial visit $k$ ($k=0,1,\ldots$) in trial $m$
($m=0,\ldots,{n_{v}}-1$). Also, we shall use the overbar to denote variable
history, for example, $\overline{A}_{m,k}=(A_{m,0},\ldots,A_{m,k})$ denotes
the patient’s history of treatment up to trial visit $k$ in trial $m$, and
define $\overline{0}=(0,\dots,0)$ and $\overline{1}=(1,\ldots,1)$. We assume
the temporal ordering of $L_{m,k},A_{m,k},Y_{m,k},C_{m,k}$ within
$[t_{m+k},t_{m+k+1})$, $\forall m,k$.

To describe the causal estimands in the next section, we shall define
potential variables $Y_{m,k}^{\overline{a}_{k}}$ and
$L_{m,k}^{\overline{a}_{k}}$. For a patient enrolled in trial $m$ (i.e., with
$E_{m}=1$) and for ${\overline{a}_{k}}$ equal to either $\overline{0}$ or
$\overline{1}$, the potential variable $Y_{m,k}^{\overline{a}_{k}}$ is an
indicator of whether that patient would have experienced the event of interest
during the time interval $[t_{m},t_{m+k+1})$ if he/she/they had, possibly
contrary to fact, received treatment $a\in\\{0,1\\}$ since the baseline of
trial $m$, i.e., from time $t_{m}$ up to time $t_{m+k}$. Analogously,
$L_{m,k}^{\overline{a}_{k}}$ denotes the time-varying covariates the patient
would have if he/she/they had received this treatment since $t_{m}$ up to time
$t_{m+k}$. Note that $Y_{m,k}^{\overline{a}_{k}}$ and
$L_{m,k}^{\overline{a}_{k}}$ are not defined for patients ineligible for trial
$m$ with $E_{m}=0$. Thus we omit the explicit conditioning on $E_{m}=1$ when
describing the causal estimand.

###  3.2 Causal estimand and assumptions

We define the per-protocol effect in trial $m$ using the MRD as follows. The
MRD at trial visit $k$ in trial $m$ is the difference between the marginal
cumulative incidence at time $t_{m+k}$ if all patients in the population
eligible for trial $m$ were always treated by time $t_{m+k}$ and the marginal
cumulative incidence if they were not treated at all by time $t_{m+k}$
keogh_causal_2023. That is,

$\mbox{MRD}_{m}(k)=\Pr(Y_{m,k}^{\overline{a}_{k}=\overline{1}}=1)-\Pr(Y_{m,k}^{\overline{a}_{k}=\overline{0}}=1)$
(1)

Identification of (1) requires key causal assumptions: no interference,
consistency, positivity and no unmeasured confounding of treatment assignment
at the trial baseline, positivity and sequential ignorability of treatment
adherence (keogh_causal_2023). If censoring due to loss of follow-up is not
independent of the potential covariates and outcome processes, then positivity
and sequential ignorability of censoring must also hold. Details of these
assumptions are provided in Section B of the Supplementary Material.

The MRD in (1) can be written equivalently as

$\displaystyle\mbox{MRD}_{m}(k)$ $\displaystyle=$
$\displaystyle\Pr(Y_{m,k}^{\overline{a}_{k}=\overline{0}}=0)-\Pr(Y_{m,k}^{\overline{a}_{k}=\overline{1}}=0)$
(2) $\displaystyle=$ $\displaystyle{\rm
E}_{V,L_{m,0}}\left\\{\Pr(Y_{m,k}^{\overline{a}_{k}=\overline{0}}=0\mid
V,L_{m,0})-\Pr(Y_{m,k}^{\overline{a}_{k}=\overline{1}}=0\mid
V,L_{m,0})\right\\}$

The counterfactual survival probabilities in (2) can be written in terms of
counterfactual discrete-time hazards:

$\Pr(Y^{\overline{a}_{k}}_{m,k}=0)={\rm
E}_{V,L_{m,0}}\left[\left\\{1-\mbox{Pr}\left(Y_{m,0}^{a_{0}}=1\mid
V,L_{m,0}\right)\right\\}\prod_{j=1}^{k}\left\\{1-\mbox{Pr}\left(Y_{m,j}^{\overline{a}_{j}}=1\mid
Y_{m,j-1}^{\overline{a}_{j-1}}=0,V,L_{m,0}\right)\right\\}\right].$

###  3.3 Marginal Structural Models with Inverse Probability Weighting

We specify an MSM in terms of the counterfactual discrete-time hazard given
baseline covariates and describe how to use Inverse Probability of Treatment
Weights(IPTW) and Inverse Probability of Censoring Weights (IPCW) to fit this
MSM. The MSM takes the following form

$\displaystyle
g\left\\{\mbox{Pr}\left(Y_{m,k}^{\overline{a}_{k}=\overline{a}}=1\mid{Y}_{m,k-1}^{\overline{a}_{k-1}=\overline{a}}=0,V,L_{m,0}\right)\right\\}=\mu(k,m,a,V,L_{m,0};\bm{\beta}),$
(3) $\displaystyle a\in\\{0,1\\}$

The MSM in (3) could be specified and fitted separately for each emulated
trial. However, a combined analysis of the sequential trials can be more
efficient and is sometimes necessary when the number of treated patients for
some trials are small. This involves making modelling assumptions about how
the MSM parameters vary across trials. For example, Danaei et al. (2013)
allowed a trial-specific intercept term but assumed that the coefficient for
treatment was the same in all trials (danaei_observational_2013).
Specifically, we can specify a pooled logistic model

$\displaystyle\mu(k,m,a,L_{m,0},V;\bm{\beta})$
$\displaystyle=\text{logit}\left\\{\mbox{Pr}\left(Y_{m,k}^{\overline{a}_{k}=\overline{a}}=1\mid{Y}_{m,k-1}^{\overline{a}_{k-1}=\overline{a}}=0,V,L_{m,0}\right)\right\\}$
$\displaystyle=\beta_{0}(m)+\beta_{1}(k)+\beta_{2}\cdot
a+\bm{\beta}_{3}^{\text{T}}V+\bm{\beta}_{4}^{\text{T}}L_{m,0},$ (4)
$\displaystyle a\in\\{0,1\\}$

with the logit link function $g(\cdot)=\text{logit}(\cdot)$, where
$\beta_{0}(m)$ is a trial-specific intercept and $\beta_{1}(k)$ is the
baseline hazard at trial visit $k$. Parametric forms or splines can be used
for specifying $\beta_{0}(m)$ and $\beta_{1}(k)$. The coefficient of treatment
$\beta_{2}$ is assumed to be common across trials, thus allowing for borrowing
strength across trials for treatment effect estimation. In (3.3) it is also
assumed that the baseline hazard and the coefficients of baseline covariates
do not vary across trials. Interactions between trials, follow-up trial
visits, treatment and baseline covariates can be included as well. Keogh et
al. (2023) (keogh_causal_2023) discussed the MSM specification in sequential
trial emulation and recommended formal tests to be performed for the
interactions.

Given the baseline covariates $V$ and $L_{m,0}$, we assume no unmeasured
confounding of baseline treatment assignment. Thus if all patients adhere to
the treatment assigned and no censoring occurs, the MSM parameters can be
estimated by fitting the pooled logistic model in (3.3) to the observed data
in the emulated trials. In practice, not all eligible patients for trial $m$
would adhere to the treatment received at baseline. Following
hernan_observational_2008 and gran_sequential_2010, we artificially censor a
patient’s follow-up in trial $m$ at the time at which he/she/they ceased to
adhere to the treatment $A_{m,0}$ received at the start of trial $m$. This
artificial censoring can introduce selection bias if it depends on time-
varying covariates measured after the start of the trial that are also
associated with the outcome. Censoring due to loss to follow-up can likewise
induce selection bias. As mentioned in Section 1, we shall use IPW to address
these two sources of selection bias.

For addressing artificial censoring due to treatment non-adherence, a
patient’s stabilised inverse probability of treatment weight
danaei_observational_2013; daniel_methods_2013; murray_causal_2021;
robins_marginal_2000 at trial visit $k$ ($k>0$) in trial $m$ is defined as

$\displaystyle
sw^{A}_{m,k}=\frac{\prod_{j=1}^{k}\Pr(A_{m,j}=a\mid\overline{A}_{m,j-1}=\overline{a},V,L_{m,0},E_{m}=1,Y_{m,j-1}=0,C_{m,j-1}=0)}{\prod_{j=1}^{k}\Pr(A_{m,j}=a\mid\overline{A}_{m,j-1}=\overline{a},V,\overline{L}_{m,j},E_{m}=1,Y_{m,j-1}=0,C_{m,j-1}=0)},$
(5) $\displaystyle a\in\\{0,1\\}$

Here,
$\Pr(A_{m,j}=a\mid\overline{A}_{m,j-1}=\overline{a},V,\overline{L}_{m,j},E_{m}=1,Y_{m,j-1}=0,C_{m,j-1}=0)$
is the conditional probability that the patient’s treatment received at trial
visit $j$ in trial $m$ remained the same as their treatment received up to
trial visit $j-1$, conditional on their observed variables up to trial visit
$j$. The stabilisation using the numerator term of (5) aims to overcome the
issue of some patients having a probability of being adherent close to 0,
which causes their unstabilised weight (numerator equal to 1) to be
excessively large robins_marginal_2000. The baseline covariates $V$ and
$L_{m,0}$ in the numerator of (5) must also be included as covariates in the
MSM of equation (3) daniel_methods_2013.

To address censoring due to loss to follow-up, the stabilised inverse
probability of censoring weight danaei_observational_2013;
daniel_methods_2013; murray_causal_2021; robins_marginal_2000 for a patient at
trial visit $k$ ($k>0$) in trial $m$ is

$\displaystyle sw^{C}_{m,k}=\frac{\prod_{j=0}^{k-1}\Pr(C_{m,j}=0\mid
C_{m,j-1}=0,\overline{A}_{m,j}=\overline{a},V,L_{m,0},E_{m}=1,Y_{m,j}=0)}{\Pi_{j=0}^{k-1}\Pr(C_{m,j}=0\mid
C_{m,j-1}=0,\overline{A}_{m,j}=\overline{a},V,\overline{L}_{m,j},E_{m}=1,Y_{m,j}=0)},$
(6) $\displaystyle a\in\\{0,1\\}$

Here,
$\Pr(C_{m,j}=0\mid\overline{C}_{m,j-1}=0,\overline{A}_{m,j}=\overline{a},V,\overline{L}_{m,j},E_{m}=1,Y_{m,j}=0)$
is the conditional probability that the patient remains in the trial during
$[t_{m,j+1},t_{m,j+2})$ given that they had been in the trial in
$[0,t_{m,j+1})$ and their observed covariate history up to time $t_{m,j}$.
Again, whenever we condition on certain variables in the numerator of (6),
these variables must also be included as covariates in the MSM of equation
(3).

Each patient’s stabilised inverse probability of treatment and censoring
weight (SIPTCW) at trial visit $k$ in trial $m$ is

$\displaystyle sw_{m,k}^{AC}=sw_{m,k}^{A}\times sw_{m,k}^{C}.$ (7)

The models used to obtain the SIPTCW could be fitted separately in each trial
or fitted on the original observational dataset. As discussed, the number of
treated patients in a single trial might be small. Similarly, the number of
patients lost to follow-up could be small for some trials. In these scenarios,
estimating the SIPTCW by fitting the models on the original observational
dataset which includes their whole treatment and censoring history becomes
desirable. For example, Danaei et al. (2013) fitted logistic models to the
treatment and censoring data from the original observational cohort to
estimate the conditional probabilities in the numerator and denominator terms
of (6) and (5) danaei_observational_2013. Note that when fitting the SIPTCW
models, we use the observed treatment and censoring data of each patient from
the visits that correspond to the baselines of the eligible trials until the
trial visits where the patient stopped adhering to the assigned treatments or
the last trial visits. If patients are eligible for multiple trials,
duplicates of the treatment and censoring data within patients are discarded.
In this paper, we follow the approach of Danaei et al. (2013) for estimating
the SIPTCW.

Under the assumptions made in Section 3.2 and that the MSM (3) and the models
for estimating the SIPTCW are correctly specified, we can then pool the
observed data from the ${n_{v}}$ trials and fit the MSM (3) to the pooled,
artificially censored and weighted data to obtain the maximum likelihood
estinate $\bm{\hat{\beta}}$ of the MSM parameters $\bm{\beta}$ in (3).

###  3.4 Estimating the causal estimand

The MRD at trial visit $k$ in trial $m$ can be estimated using the parameter
estimates $\hat{\bm{\beta}}$ from the MSM by the empirical standardization
formula:

$\displaystyle\widehat{\mbox{MRD}}_{m}(k)$ $\displaystyle=$
$\displaystyle\frac{1}{n_{m}}\sum_{i=1}^{n}E_{m,i}\prod_{j=0}^{k}\left\\{1-g^{-1}\left\\{\mu(j,m,a=0,V_{i},L_{m0,i};\hat{\bm{\beta}})\right\\}\right\\}$
(8)
$\displaystyle-\frac{1}{n_{m}}\sum_{i=1}^{n}E_{m,i}\prod_{j=0}^{k}\left\\{1-g^{-1}\\{\mu(j,m,a=1,V_{i},L_{m0,i};\hat{\bm{\beta}})\\}\right\\}$

Here, $i$ indexes the patient in the original observational study, and
$n_{m}=\sum_{i=1}^{n}E_{m,i}$ is the total number of patients enrolled in the
$m$th trial. $\widehat{\mbox{MRD}}_{m}(k)$ is an estimate of the MRD in the
population of patients eligible for trial $m$. Alternatively, we could
standardise to a population with a different distribution of $(V,L_{m0})$.

To summarise, the MRD in (2) is estimated by the following steps:

  1. 1.

Estimate the inverse probability of treatment and censoring weights using data
from the observational cohort.

  2. 2.

Expand the observational data by assigning patients to eligible sequential
trials with corresponding follow-up trial visits and artificially censor
patients’ follow-up when they were no longer adhering to the treatment
assigned at the baseline of each emulated trial.

  3. 3.

Estimate the MSM for the sequential trials using the expanded, artificially
censored data with the estimated weights.

  4. 4.

Create two datasets containing patients from a target population with their
baseline covariates data, setting their treatment assignment to either
treatment arm, and calculate the fitted survival probability of each patient
in each of these two datasets.

  5. 5.

Estimate the MRD by averaging the survival probabilities in each of the two
datasets and taking the difference between these two averages, as is done in
equation (8).

In this paper, we use the R package TrialEmulation rezvani_trialemulation_2023
to implement these steps. More details about the TrialEmulation package can be
found in Section D of the Supplementary Material.

##  4 Constructing confidence intervals of the marginal risk difference

In this section, we describe the methods for constructing CIs for the MRD
based on the sandwich variance estimator and the nonparametric bootstrap, the
proposed method based on the LEF bootstrap and Jackknife resampling.

###  4.1 Sandwich variance estimator

The sandwich variance estimator accounts for the inverse probability weights
and the correlation induced by patients being eligible to multiple trials
(hernan_marginal_2000; danaei_observational_2013). Specifically, for the
weighted pooled logistic model such as in (3.3), the sandwich variance
estimator of $\hat{\bm{\beta}}$ is

$\hat{\bm{\Sigma}}=\left\\{\sum_{i=1}^{n}\frac{\partial{\bf
U}_{i}(\bm{\beta})}{\partial\bm{\beta}^{\text{T}}}\right\\}_{\bm{\beta}=\hat{\bm{\beta}}}^{-1}\left\\{\sum_{i=1}^{n}{\bf
U}_{i}(\hat{\bm{\beta}}){\bf
U}_{i}(\hat{\bm{\beta}})^{\text{T}}\right\\}\left\\{\frac{\partial{\bf
U}_{i}(\bm{\beta})}{\partial\bm{\beta}^{\text{T}}}\right\\}_{\bm{\beta}=\hat{\bm{\beta}}}^{-1}$

where ${\bf U}_{i}(\hat{\bm{\beta}})$ is the weighted score function of the
pooled logistic model evaluated at $\hat{\bm{\beta}}$ for patient $i$.

We follow the parametric bootstrap algorithm of Mandel (2013) to construct
simulation-based CIs of the MRD as follows mandel_simulation-based_2013:

  1. 1)

Obtain the parameter estimate $\hat{\bm{\beta}}$ and the sandwich variance
matrix estimate $\hat{\bm{\Sigma}}$ of the MSM in (3.3).

  2. 2)

Draw an i.i.d. sample $\bm{\beta}^{(1)},...,\bm{\beta}^{(S)}$ of size $S$ (say
$S=1000$) from the multivariate Normal distribution
$\mbox{MVN}(\hat{\bm{\beta}},\hat{\bm{\Sigma}})$ with mean $\hat{\bm{\beta}}$
and variance-covariance $\hat{\bm{\Sigma}}$ matrix.

  3. 3)

For each vector $\bm{\beta}^{(s)}$ ($s=1,...,S$), estimate the MRD at each
trial visit by setting $\bm{\beta}^{(s)}$ as the MSM parameters.

  4. 4)

Use the 2.5th and 97.5th percentiles of these $S$ MRD estimates at each trial
visit as the lower and upper-bounds of the 95% CI.

This is the only CI method implemented in the TrialEmulation package.

###  4.2 Nonparametric bootstrap with the pivot method

We used the non-Studentized pivot method carpenter_bootstrap_2000 to construct
CIs based on nonparametric bootstrap. The rationale behind the non-Studentized
pivot CI is that the relationship between the true MRD and the MRD estimate
based on the original data is approximately the same as the relationship
between the MRD estimate based on the original data and the MRD estimate based
on a bootstrap sample banks_lecture_2017. Specifically, we followed the steps
described below:

  1. 1)

Draw $B$ bootstrap samples from the observational cohort data, treating the
$n$ patients as the resampling units.

  2. 2)

For each bootstrap sample $b$, $b=1,\ldots,B$ obtain the bootstrap parameter
estimate $\hat{\bm{\beta}}^{(b)}$ and estimate the MRD at each trial visit $k$
in trial $m$ ${\mbox{MRD}}_{m}(k)$ using the method in Section 3.

  3. 3)

Define the lower and upper bounds of the 95% CI for the MRD at each trial
visit $k$ in trial $m$ as, respectively,
$2\widehat{\mbox{MRD}}_{m}(k)-\widehat{\mbox{MRD}}_{m}(k)^{*}_{(0.975)}$ and
$2\widehat{\mbox{MRD}}_{m}(k)-\widehat{\mbox{MRD}}_{m}(k)^{*}_{(0.025)},$
where $\widehat{\mbox{MRD}}_{m}(k)$ is the point estimate of the MRD at trial
visit $k$ in trial $m$ estimated from the original data, and
$\widehat{\mbox{MRD}}_{m}(k)^{*}_{(0.025)}$,
$\widehat{\mbox{MRD}}_{m}(k)^{*}_{(0.975)}$ are the 2.5th and 97.5th
percentiles of the $B$ bootstrap MRD estimates at trial visit $k$ in trial
$m$, respectively.

###  4.3 Linearised Estimating Function bootstrap

Linearised Estimating Function (LEF) bootstrap is a computationally efficient
bootstrap method proposed by Binder et al. (2004) binder_design-based_2004. It
is an adaptation of the estimating function bootstrap developed by Hu and
Kalbfleisch hu_estimating_2000. The main advantages of LEF bootstrap over
nonparametric bootstrap are reduced computational time and convergence issues.
In terms of computational time, unlike the nonparametric bootstrap, which
involves fitting a regression model to each bootstrap sample, the LEF
bootstrap method requires fitting this model only once: to the original
dataset. This can reduce computational time considerably, especially when
iterative procedures are used to fit the model. In terms of convergence
issues, Binder et al. (2004) found that when using nonparametric bootstrap for
logistic regression it was possible to have several bootstrap samples for
which the parameter estimation algorithm would not converge, due to ill-
conditioned matrices that were not invertible binder_design-based_2004. For
both these reasons, the LEF bootstrap method may have advantages in our
sequential trial emulation setting, where logistic regression models are used
both for estimating the SIPTCWs and for fitting MSMs. We now explain how the
LEF bootstrap works in the general situation where the goal is to construct a
CI for some function of a parameter vector $\theta$. In Section 4.3.1, we
shall describe how to extend this general method to the specific setting of
sequential trial emulation.

Let ${\bf U}(\mbox{\boldmath{$\theta$}})$ denote the estimating function for
$\theta$ (note that ${\bf U}$ here is different from the ${\bf U}_{i}$ in the
previous section). Let ${\bf U}^{\rm org}(\mbox{\boldmath{$\theta$}})$ denote
the sum of ${\bf U}(\mbox{\boldmath{$\theta$}})$ over the $n$ patients in the
original dataset. Then ${\bf U}^{\rm
org}(\mbox{\boldmath{$\theta$}})=\mathbf{0}$ is the estimating equation for
$\theta$ based on the original dataset. Let $\hat{\mbox{\boldmath{$\theta$}}}$
denote the estimate of $\theta$ obtained by solving this equation. Now suppose
that $B$ bootstrap samples have been generated by resampling with replacement
from this original dataset. Let ${\bf
U}^{(b)}(\mbox{\boldmath{$\theta$}})=\mathbf{0}$ denote the estimating
equation for $\theta$ based on the $b$th bootstrap dataset, and let
$\hat{\mbox{\boldmath{$\theta$}}}^{(b)}$ denote the corresponding estimate of
$\theta$ ($b=1,\ldots,B$). If we (Taylor) linearise the function ${\bf
U}^{(b)}(\mbox{\boldmath{$\theta$}})$ around
$\hat{\mbox{\boldmath{$\theta$}}}$, we obtain

${\bf U}^{(b)}({\mbox{\boldmath{$\theta$}}})\approx{\bf
U}^{(b)}(\hat{\mbox{\boldmath{$\theta$}}})+\left\\{\frac{\partial{\bf
U}^{(b)}({\mbox{\boldmath{$\theta$}}})}{\partial\mbox{\boldmath{$\theta$}}}\right\\}_{\mbox{\boldmath{$\theta$}}=\hat{\mbox{\boldmath{$\theta$}}}}(\mbox{\boldmath{$\theta$}}-\hat{\mbox{\boldmath{$\theta$}}}).$

From this and the fact that ${\bf
U}^{(b)}(\hat{\mbox{\boldmath{$\theta$}}}^{(b)})=\mathbf{0}$ by definition, we
obtain

$\mathbf{0}\approx{\bf
U}^{(b)}(\hat{\mbox{\boldmath{$\theta$}}})+\left\\{\frac{\partial{\bf
U}^{(b)}({\mbox{\boldmath{$\theta$}}})}{\partial\mbox{\boldmath{$\theta$}}}\right\\}_{\mbox{\boldmath{$\theta$}}=\hat{\mbox{\boldmath{$\theta$}}}}(\hat{\mbox{\boldmath{$\theta$}}}^{(b)}-\hat{\mbox{\boldmath{$\theta$}}})$

Rearranging the terms, we get

$\hat{\mbox{\boldmath{$\theta$}}}^{(b)}\approx\hat{\mbox{\boldmath{$\theta$}}}-\left\\{\frac{\partial{\bf
U}^{(b)}({\mbox{\boldmath{$\theta$}}})}{\partial\mbox{\boldmath{$\theta$}}}\right\\}_{\mbox{\boldmath{$\theta$}}=\hat{\mbox{\boldmath{$\theta$}}}}^{-1}{\bf
U}^{(b)}(\hat{\mbox{\boldmath{$\theta$}}}).$

Replacing the matrix $\left\\{\partial{\bf
U}^{(b)}(\mbox{\boldmath{$\theta$}})/\partial\mbox{\boldmath{$\theta$}}\right\\}_{\mbox{\boldmath{$\theta$}}=\hat{\mbox{\boldmath{$\theta$}}}}^{-1}$
by $\left\\{\partial{\bf
U}(\mbox{\boldmath{$\theta$}})/\partial\mbox{\boldmath{$\theta$}}\right\\}_{\mbox{\boldmath{$\theta$}}=\hat{\mbox{\boldmath{$\theta$}}}}^{-1}$,
we obtain the following approximation of $\hat{\bm{\theta}}^{(b)}$:

$\displaystyle\hat{\mbox{\boldmath{$\theta$}}}^{(b)}_{\mbox{LEF}}$
$\displaystyle\approx$
$\displaystyle\hat{\mbox{\boldmath{$\theta$}}}-\left\\{\frac{\partial{\bf
U}({\mbox{\boldmath{$\theta$}}})}{\partial\mbox{\boldmath{$\theta$}}}\right\\}_{\mbox{\boldmath{$\theta$}}=\hat{\mbox{\boldmath{$\theta$}}}}^{-1}{\bf
U}^{(b)}(\hat{\mbox{\boldmath{$\theta$}}})=\hat{\mbox{\boldmath{$\theta$}}}+\text{vcov}({\hat{\bm{\theta}}}){\bf
U}^{(b)}(\hat{\mbox{\boldmath{$\theta$}}}),$ (9)

where $\text{vcov}({\hat{\bm{\theta}}})$ is the model-based variance matrix
based on the original dataset. A CI for any given function of $\theta$ can now
be constructed in the same way as is done in Section 4.2 for the nonparametric
bootstrap, except with $\hat{\mbox{\boldmath{$\theta$}}}^{(b)}_{\mbox{LEF}}$
used in place of $\hat{\bm{\beta}}^{(b)}$. This approximation of
$\hat{\bm{\beta}}^{(b)}$ by
$\hat{\mbox{\boldmath{$\theta$}}}^{(b)}_{\mbox{LEF}}$ avoids the need to solve
the estimating equations ${\bf U}^{(b)}(\mbox{\boldmath{$\theta$}})=0$ for
each $b=1,\ldots,B$ (as would be done if applying nonparametric bootstrap).

We propose the following two ways of using the LEF bootstrap approach to
construct CIs for the MRD.

####  4.3.1 Approach 1: LEF bootstrap for the MSM parameters

This approach applies the Taylor linearisation only to the estimating function
of the pooled logistic regression for the MSM (3.3), with the SIPTCWs first
being estimated from each bootstrap dataset by fitting the weight models to
that bootstrap dataset, as is done in nonparametric bootstrap. We shall use
$w_{m,k,i}$ to denote the estimated SIPTCW for patient $i$ ($i=1,\ldots,n$) in
the original dataset at trial visit $k$ in trial $m$ (if patient $i$ is not
enrolled in trial $m$, i.e. $E_{m,i}=0$, then $w_{m,k,i}=0$). Analogously,
$w^{(b)}_{m,k,i}$ will denote the estimated SIPTCW for patient $i$ at trial
visit $k$ in trial $m$ in the $b$th bootstrap dataset. The procedure is as
follows.

  1. (1)

Estimate the weights $w_{m,k,i}$ using the original data. Using these
estimated weights $w_{m,k,i}$, fit the weighted pooled logistic model to the
original dataset and obtain the point estimate $\hat{\bm{\beta}}$ of the MSM
parameters.

  2. (2)

Create $B$ bootstrap samples using the $n$ patients as resampling units. Using
the $b^{\text{th}}$ bootstrap dataset, estimate the SIPTCWs $w^{(b)}_{m,k,i}$
for the patients in bootstrap dataset $b$ ($b=1,...,B$).

  3. (3)

Calculate approximate bootstrap parameter estimate $\hat{\bm{\beta}}_{\rm
LEF}^{(b)}$ from the $b$th bootstrap sample using the weights
$w^{(b)}_{m,k,i}$ according to the following formula:

$\displaystyle\hat{\bm{\beta}}_{\rm LEF}^{(b)}$
$\displaystyle=\hat{\bm{\beta}}-\left(\frac{\partial{\bf
U}(\bm{\beta})}{\partial\bm{\beta}}\right)_{\bm{\beta}=\hat{\bm{\beta}}}^{-1}\hat{\bf
U}^{(b)}(\hat{\bm{\beta}})$
$\displaystyle=\hat{\bm{\beta}}+\text{vcov}({\hat{\bm{\beta}}})\hat{\bf
U}^{(b)}(\hat{\bm{\beta}})$

For our case with a pooled logistic model, where
$g^{-1}\\{\mu_{m,k,i}({\bm{\beta}})\\}=\mbox{Pr}(Y_{m,k,i}=1\mid
Y_{m,k-1,i}=0,V_{i},L_{m,0,i},E_{m,i}=1;{\bm{\beta}})$ with
$g^{-1}(\cdot)=\text{expit}(\cdot)$, the inverse of the logit link function
$g(\cdot)$, this formula can be written as

$\displaystyle\hat{\bm{\beta}}_{\rm LEF}^{(b)}=\hat{\bm{\beta}}$
$\displaystyle-$
$\displaystyle\left(\sum_{m=0}^{{n_{v}}-1}\sum_{i=1}^{n}E_{m,i}\sum_{k=0}^{q_{m,i}-1}\;w_{m,k,i}\;g^{-1}\\{\mu_{m,k,i}(\hat{\bm{\beta}})\\}\left[1-g^{-1}\\{\mu_{m,k,i}(\hat{\bm{\beta}})\\}\right]{\bf
X}_{m,k,i}{\bf X}_{m,k,i}^{\text{T}}\right)^{-1}$ (10)
$\displaystyle\sum_{m=0}^{{n_{v}}-1}\sum_{i=1}^{n}E^{(b)}_{m,i}\sum_{k=0}^{q_{m,i}-1}w^{(b)}_{m,k,i}\;{\bf
X}^{(b)}_{m,k,i}\;\left\\{Y^{(b)}_{m,k,i}-g^{-1}\\{\mu^{(b)}_{m,k,i}(\hat{\bm{\beta}})\\}\right\\},$

where $q_{m,i}$ is the total number of trial visits made by patient $i$
eligible in trial $m$ before being artificially censored, loss to follow-up,
the occurrence of outcome event or reaching the end of trial $m$,
$g^{-1}\\{\mu_{m,k,i}(\hat{\bm{\beta}})\\}$ is the estimated discrete-time
hazard at trial visit $k$ for patient $i$ in trial $m$ evaluated at
$\hat{\bm{\beta}}$ using the MSM (3.3) fitted to the original data, ${\bf
X}_{m,k,i}$ is the corresponding design vector for the discrete-time hazard
MSM in the original data, ${\bf X}^{(b)}_{m,k,i}$ is the design vector for the
discrete-time hazard MSM at visit $k$ for patient $i$ in trial $m$ in the
$b^{\text{th}}$ bootstrap sample, $Y^{(b)}_{m,k,i}$ is the observed outcome
indicator at trial visit $k$ for patient $i$ in trial $m$ in the
$b^{\text{th}}$ bootstrap sample, and
$g^{-1}\\{\mu^{(b)}_{m,k,i}(\hat{\bm{\beta}})\\}$ is the estimated discrete-
time hazard at trial visit $k$ for patient $i$ in trial $m$ evaluated at
$\hat{\bm{\beta}}$ in the $b^{\text{th}}$ bootstrap sample. $E_{m,i}$ is the
patient’s eligibility indicator for trial $m$ in the original data, and
$E^{(b)}_{m,i}$ is the patient’s eligibility indicator for trial $m$ in the
$b^{\text{th}}$ bootstrap sample.

  4. (4)

For each bootstrap sample $b$, estimate the MRD at each trial visit by setting
$\hat{\bm{\beta}}_{\rm LEF}^{(b)}$ as the MSM parameters.

  5. (5)

Construct the pivot CI using the $2.5$th and $97.5$th percentiles of the $B$
LEF bootstrap estimates of the MRD at each trial visit, as in Section 4.2.

####  4.3.2 Approach 2: LEF bootstrap for the model parameters for SIPTCWs and
the MSM parameters

In the second approach, the Taylor linearisation is applied to the estimating
functions of both the models for estimating SIPTCWs and the pooled logistic
regression for the MSM (3.3). Approach 2 should theoretically be even more
computationally efficient than Approach 1 as no matrix inversion is required
for the models of SIPTCWs either. Moreover, Approach 2 could be potentially
useful when many covariates are included for estimating SIPTCWs using logistic
models and ill-conditioned matrices can arise for some bootstrap samples.

Because multiple models have to be fitted to obtain the SIPTCWs, we explain
the steps for implementing the LEF bootstrap using the model for the
denominator term of (5) given $\overline{A}_{m,k-1}=\overline{1}$ as an
illustration. Implementing the LEF bootstrap is analogous for all other models
involved in (5) and (6).

Let
$p_{m,k}=\Pr(A_{m,k}=1\mid\overline{A}_{m,k-1}=\overline{1},V,L_{m,k},E_{m}=1,Y_{m,k-1}=0,C_{m,k-1}=0)$,
the conditional probability that the patient remains treated at trial visit
$k$ in trial $m$ given they received treatment up to trial visit $k-1$,
conditional on their observed variables up to trial visit $k$. Suppose that a
logistic regression model is assumed for $p_{m,k}$,

$\displaystyle\text{logit}\,(p_{m,k})={\bf Z}_{m,k}^{\text{T}}{\bm{\gamma}},$
(11)

where ${\bf Z}_{m,k}$ is the design vector and ${\bm{\gamma}}$ is the
regression parameter vector. We note that the design matrix ${\bf Z}_{m,k}$
for this regression model only contains (transformations of) observations from
patients who were always treated in trial $m$ up until the previous trial
visit $k-1$, i.e. $\overline{A}_{m,k-1}=\overline{1}$ . The procedure is as
follows.

  1. (1)

Fit the weighted pooled logistic model to the original data and obtain the
point estimate $\hat{\bm{\beta}}$ of the MSM parameters, using the estimated
weights $w_{m,k,i}$ based on the original data.

  2. (2)

Create $B$ bootstrap samples with patients as resampling units, and assign a
bootstrap sampling weight $s^{(b)}_{i}$, i.e., the number of times that a
patient is sampled in the $b^{\text{th}}$ bootstrap sample ($b=1,\dots,B$),
for each patient (e.g., $s^{(b)}_{i}$ = 3 if patient $i$ was sampled 3 times
in the $b^{\text{th}}$ bootstrap sample).

  3. (3)

Calculate the LEF bootstrap parameter estimates of the models for estimating
SIPTCWs for each bootstrap sample according to (9). For example, the LEF
bootstrap estimates $\hat{\bm{\gamma}}^{(b)}$ of $\bm{\gamma}$ in (11) can be
obtained by

$\displaystyle\hat{\bm{\gamma}}^{(b)}=\hat{\bm{\gamma}}$ $\displaystyle-$
$\displaystyle\left[\sum_{m=0}^{n_{v}-1}\sum_{i=1}^{n}E_{m,i}\sum_{k=0}^{q_{m,i}-1}\mathds{1}_{\\{\overline{A}_{m,k-1,i}=\overline{1}\\}}p_{m,k,i}(\hat{\bm{\gamma}})\left\\{1-p_{m,k,i}(\hat{\bm{\gamma}})\right\\}{\bf
Z}_{m,k,i}{\bf Z}_{m,k,i}^{\text{T}}\right]^{-1}$ (12)
$\displaystyle\sum_{m=0}^{n_{v}-1}\sum_{i=1}^{n}E^{(b)}_{m,i}\sum_{k=0}^{q_{m,i}-1}\mathds{1}_{\\{\overline{A}^{(b)}_{m,k-1,i}=\overline{1}\\}}s^{(b)}_{i}{\bf
Z}^{(b)}_{m,k,i}\left\\{A^{(b)}_{m,k,i}-p^{(b)}_{m,k,i}(\hat{\bm{\gamma}})\right\\},$

where $q_{m,i}$ is the total number of trial visits made by patient $i$
eligible in trial $m$ before being artificially censored, loss to follow-up,
the occurrence of outcome event or reaching the end of trial $m$,
$\mathds{1}_{\\{\overline{A}_{m,k-1,i}=\overline{1}\\}}$ is an indicator for
whether patient $i$ has always been treated up to trial visit $k-1$ in trial
$m$, $p_{m,k,i}(\hat{\bm{\gamma}})$ is the estimated probability of treatment
adherence given previous treatment at trial visit $k$ for patient $i$ in trial
$m$ evaluated at $\hat{\bm{\gamma}}$ using the logistic regression model (11)
fitted to the original data, ${\bf Z}_{m,k,i}$ is the corresponding design
vector for the logistic regression model in the original data,
$\mathds{1}_{\\{\overline{A}^{(b)}_{m,k-1,i}=\overline{1}\\}}$ is an indicator
for whether patient $i$ has always been treated up to trial visit $k-1$ in
trial $m$ in the $b^{\text{th}}$ bootstrap sample, ${\bf Z}^{(b)}_{m,k,i}$ is
the design vector for the logistic regression model at visit $k$ for patient
$i$ in trial $m$ in the $b^{\text{th}}$ bootstrap sample, $A^{(b)}_{m,k,i}$ is
the observed treatment indicator at trial visit $k$ for patient $i$ in trial
$m$ in the $b^{\text{th}}$ bootstrap sample, and
$p^{(b)}_{m,k,i}(\hat{\bm{\gamma}})$ is the estimated probability of treatment
adherence given previous treatment at trial visit $k$ for patient $i$ in trial
$m$ evaluated at $\hat{\bm{\gamma}}$ in the $b^{\text{th}}$ bootstrap sample.
$E_{m,i}$ is the patient’s eligibility indicator for trial $m$ in the original
data, and $E^{(b)}_{m,i}$ is the patient’s eligibility indicator for trial $m$
in the $b^{\text{th}}$ bootstrap sample.

  4. (4)

Based on LEF bootstrap parameter estimates for the models for estimating
SIPTCWs, we calculate a new set of weights $\tilde{w}^{(b)}_{m,k,i}$ for the
$b^{\text{th}}$ bootstrap sample.

  5. (5)

Replace ${w}^{(b)}_{m,k,i}$ with $\tilde{w}^{(b)}_{m,k,i}$ in (10) and obtain
the $b^{\text{th}}$ bootstrap parameter estimate $\hat{\bm{\beta}}_{\rm
LEF}^{(b)}$ using (10) .

  6. (6)

For each bootstrap sample $b$, estimate the MRD at each trial visit by setting
$\hat{\bm{\beta}}^{(b)}$ as the MSM parameters.

  7. (7)

Construct the pivot CI using the $2.5$th and $97.5$th percentiles of the $B$
LEF bootstrap estimates of the MRD at each trial visit as in Section 4.2.

Steps (5)-(7) are exactly the same as Steps (3)-(5) in Approach 1, except with
different sets of weights $\tilde{w}^{(b)}_{m,k,i}$.

###  4.4 Jackknife

We used Jackknife resampling to construct two types of Jackknife CIs: 1) using
a Jackknife estimate of the MRD standard error tibshirani_introduction_1994,
or 2) using the Jackknife variance estimator lipsitz_using_1990;
friedl_jackknife_2001 of $\hat{\bm{\beta}}$. For longitudinal studies, the
sampling unit is the patient ID and for from a dataset with $n$ patients
observed, $n$ Jackknife samples are obtained by removing all of the $i$th
patient’s observations lipsitz_using_1990; friedl_jackknife_2001 in each
sample ($i=1,\dots,n$). Adapting this to sequential trial emulation, this
translates to removing all of the $i$th patient’s observations from all trials
they are eligible to in each Jackknife sample.

####  4.4.1 Approach 1: Wald-type confidence interval using the Jackknife
estimate of MRD standard error

We use the Jackknife standard error estimate of the MRD and a Normal
approximation of the distribution of the MRD to obtain CIs as follows:

  1. 1)

Obtain the parameter estimate $\hat{\bm{\beta}}$ of the MSM in (3.3) and the
MRD estimate $\widehat{\mbox{MRD}}_{m}(k)$ at each trial visit $k$ in trial
$m$ using the method in Section 3.

  2. 2)

Draw $n$ leave-one-out Jackknife samples from the observational cohort data,
treating the $n$ patients as the resampling units.

  3. 3)

For each of the $n$ Jackknife samples where we leave out patient $i$’s
observations, obtain the Jackknife parameter estimate
$\hat{\bm{\beta}}^{(-i)}$ and estimate $\mbox{MRD}^{(-i)}_{m}(k)$ at each
trial visit $k$ in trial $m$ using the method in Section 3.

  4. 4)

Obtain the Jackknife standard error estimate $\widehat{SE}^{J}_{m}(k)$ of the
MRD at trial visit $k$ in trial $m$,

$\displaystyle\widehat{SE}^{J}_{m}(k)=\left[\frac{n-1}{n}\sum_{i=1}^{n}\left(\widehat{\mbox{MRD}}^{(-i)}_{m}(k)-\widehat{\mbox{MRD}}_{m}(k)\right)^{2}\right]^{1/2}$
(13)

where $\widehat{\mbox{MRD}}_{m}(k)$ is the point estimate of the MRD at trial
visit $k$ in trial $m$ estimated from the original data.

  5. 5)

Define the lower and upper bounds of the 95% CI for the MRD at trial visit $k$
in trial $m$, respectively,
$\widehat{\mbox{MRD}}_{m}(k)-1.96\cdot\widehat{SE}^{J}_{m}(k)$ and
$\widehat{\mbox{MRD}}_{m}(k)+1.96\cdot\widehat{SE}^{J}_{m}(k).$

####  4.4.2 Approach 2: Multivariate Normal sampling using the Jackknife
Variance Estimator of MSM parameters

We use the Jackknife variance estimator lipsitz_using_1990;
friedl_jackknife_2001 of $\hat{\bm{\beta}}$ to construct simulation-based CIs
of the MRD in a manner similar to the sandwich variance estimator-based CIs as
follows:

  1. 1)

Obtain the parameter estimate $\hat{\bm{\beta}}$ of the MSM in (3.3).

  2. 2)

Draw $n$ leave-one-out Jackknife samples from the observational cohort data,
treating the $n$ patients as the resampling units.

  3. 3)

For each of the $n$ Jackknife samples where we leave out patient $i$’s
observations, obtain the Jackknife parameter estimate
$\hat{\bm{\beta}}^{(-i)}$.

  4. 4)

Obtain the Jackknife variance matrix estimate $\hat{\bm{V_{J}}}$ of the MSM:

$\displaystyle\hat{\bm{V_{J}}}=\frac{1}{n(n-1)}\sum_{i=1}^{n}(b_{i}-\bar{\beta})(b_{i}-\bar{\beta})^{T}$
(14)

where $b_{i}=n\hat{\bm{\beta}}-(n-1)\hat{\bm{\beta}}^{(-i)}$ for $i=1,...,n$
and $\bar{\beta}=n^{-1}\sum_{i=1}^{n}b_{i}$.

  5. 5)

Construct a 95% CI by repeating Steps (2)-(4) of Section 4.1, replacing the
sandwich variance matrix estimate $\hat{\bm{\Sigma}}$ with $\hat{\bm{V_{J}}}$.

##  5 Simulation study

We conducted a simulation study to compare the performance of 95% CIs for the
MRD in per-protocol analysis of sequentially emulated trials using the
sandwich variance estimator, nonparametric bootstrap, the proposed LEF
bootstrap and Jackknife methods. For better quantification of the influence of
different setups for the treatment and outcome processes, we focused on
scenarios without dependent censoring due to loss to follow-up.

###  5.1 Study setup

####  5.1.1 Data generating mechanism

We used the algorithm described in Young and Tchetgen Tchetgen (2014)
young_simulation_2014 to simulate data with interrelationships between time-
varying confounders, treatment assignments and the survival time that arise
typically in an observational cohort. Specifically, previous treatments would
affect time-varying confounders that are also associated with current
treatments received and the risk of the survival outcome. The data generating
mechanism is described in Table 2.

Table 2: Summary of data generating mechanism of the simulation study. Data simulation setting specifications |  $n$: number of patients
---|---
|  $n_{v}=5$, number of visits
|  $t_{j}=0,...,n_{v}-1$: visit time for visit $j$
|  $\alpha_{a}$: intercept in the treatment model, representing the baseline
rate of treatment initiation
|  $\alpha_{c}$: coefficient that describes the strength of confounding due to
time-varying variable $X_{1,t_{j}}$ on treatment assignment and outcome
|  $\alpha_{y}$: intercept term in the discrete time-hazard model,
representing the baseline hazard
Time-varying confounder |  $X_{1,t_{j}}\sim N(Z_{t_{j}}-0.3A_{t_{j-1}},1)$, with $A_{t_{-1}}=0$.
Time-invariant confounder |  $X_{2}\sim N(0,1)$
Treatment |  $\text{logit}\\{\Pr(A_{t_{j}}=1\mid A_{t_{j-1}},X_{1,t_{j}},X_{2},Y_{t_{j-1}}=0)\\}$
|  $=\alpha_{a}+0.05A_{t_{j-1}}+\alpha_{c}X_{1,t_{j}}+0.2X_{2}$, where
$Y_{t_{-1}}=0$.
Discrete-time hazard of the survival outcome |  $\text{logit}\\{\Pr(Y_{t_{j}}=1\mid A_{t_{j}},X_{1,t_{j}},X_{2},Y_{t_{j-1}}=0)\\}$
|  $=\alpha_{y}-0.5A_{t_{j}}+\alpha_{c}X_{1,t_{j}}+X_{2}$
Censoring |  No censoring
Trial Eligibility |  $E_{t_{j}}=1$ if patient has not received treatment before $t_{j}$
|  and has not experienced the outcome event before $t_{j}$;
|  $E_{t_{j}}=0$ otherwise

####  5.1.2 Monte Carlo simulation settings

We considered three settings with different outcome event rates, by setting
the baseline hazards $\alpha_{y}$ as $-4.7$, $-3.8$ and $-3$ to emulate low
($5-6.5\%$), moderate ($10-14\%$) and high ($20-25\%$) percentage of patients
experiencing the event during follow-up in the simulated data, respectively.
In total, we investigated 81 different scenarios for generating the simulated
data using the mechanism in Table 2, by considering the combinations of the
specifications presented in Table 3. We varied the number of patients,
confounding strength of a time-varying confounder and treatment prevalence
(i.e., percentage of patients who ever received treatments during follow-up).
By varying the intercept term of the treatment model $\alpha_{a}$, we could
generate a low ($25-30\%$), medium ($50-60\%$) and high ($75-80\%$) treatment
prevalence in the simulated data. For each scenario, we generated 1000
simulated datasets.

Table 3: Summary of the specifications for the 81 scenarios considered in the simulations. Outcome even rate | Sample size | Confounding strength | Treatment prevalence
---|---|---|---
Low: $\alpha_{y}=-4.7$ | Small: $n=200$ | $\alpha_{c}=0.1$ | Low: $\alpha_{a}=-1$
Medium: $\alpha_{y}=-3.8$ | Medium: $n=1000$ | $\alpha_{c}=0.5$ | Medium: $\alpha_{a}=0$
High: $\alpha_{y}=-3$ | Large: $n=5000$ | $\alpha_{c}=0.9$ | High: $\alpha_{a}=1$

An example of the number of events we can observe following the data
generating algorithm in each treatment arm at each trial visit for each data
scenario can be found in Tables 10, 11 and 12 of the Supplementary Materials.

####  5.1.3 Estimation and inference

For each simulated dataset, we emulated 5 trials ($m=0,\ldots,4$), with trial
0 including trial visits $k=0,\ldots,4$, trial 1 including trial visits
$k=0,\ldots,3$, and so on. Our estimand of interest is the MRDs at trial
visits $k=0,\ldots,4$ for patients eligible in trial $m=0$. We chose trial 0
patients as the target population because they had the longest length of
follow-up so we can assess the methods for estimation and inference of the
MRDs at later visits.

Correctly specified logistic regression models were fitted to the simulated
data to estimate the denominator terms of the stabilised IPTWs (5), while we
estimated the numerator terms by fitting an intercept-only logistic model
stratified by treatments received at the immediately previous visit
$A_{t_{j-1}}$. Then the following weighted pooled logistic model was fitted to
the artificially censored data in the five sequential trials to estimate the
counterfactual discrete-time hazard:

$\displaystyle\mu(k,m,a,L_{m,0},V;\bm{\beta})$
$\displaystyle=\text{logit}\\{\Pr(Y_{m,k}=1\mid
Y_{m,k-1}=0,X_{1m,0},X_{2})\\}$ $\displaystyle=\beta_{0,k}+\beta_{1,k}\cdot
a+\beta_{2,k}X_{1,m,0}+\beta_{3,k}X_{2}$

where $X_{1,m,0}$ is the value of the time-varying confounder at the baseline
of trial $m$ ($m=0,\ldots,4$) and $\beta_{0,k}$, $\beta_{1,k}$, $\beta_{2,k}$
$\beta_{3,k}$ are regression coefficients that vary by trial visit number
($k=0,\ldots,4$) but are assumed to be the same across trials. Note that we
were not able to correctly specify the MSM for counterfactual discrete-time
hazard via pooled logistic regression models, due to the non-collapsibility of
the logistic model used in the data-generating mechanism (keogh_causal_2023).
To reduce potential bias due to this misspecification, we specified a rich
model by allowing the coefficients of treatment and confounders to vary by
trial visits. Based on Model (5.1.3), we used the methods in Section 3.4 to
obtain the point estimates of the MRD for the population of patients eligible
for trial $0$.

Finally, we applied the methods for constructing 95% CIs of the MRD in Section
4 to each simulated dataset, where we set $S=500$ for the CIs based on the
sandwich variance estimator and $B=500$ for the CIs based on the bootstrap
methods. We restricted our simulation studies to applying the two Jackknife CI
construction methods to simulated datasets with $n=200$ only, as it would be
computationally inefficient to use the Jackknife CI construction methods for
larger sample sizes. A pseudo-code algorithm for the simulation study is
presented in Section C of the Supplementary Materials.

####  5.1.4 True values

True values of the MRDs for each simulation scenario were obtained, as
proposed by Keogh et al. keogh_causal_2023, by generating data for a very
large randomized controlled trial. The true marginal risks of patients in
trial 0 were approximated by Kaplan-Meier estimates from two extremely large
datasets with 1,000,000 patients, which were generated by enforcing the two
treatment strategies, the ‘always treated’ (by setting $A_{t_{j}}=1$ for
$j=0,\ldots,4$) and ‘never treated’ (by setting $A_{t_{j}}=0$ for
$j=0,\ldots,4$) in the same data generation algorithm in Table 2. The true MRD
was obtained by taking the difference between these counterfactual marginal
risks.

####  5.1.5 Performance measures

Empirical coverage of the CIs based on different methods was calculated for
each trial visit in each simulation scenario. This was done by dividing the
number of times that a CI of the MRD contained the true MRD by 1000 minus the
number of times there was an error in the CI construction. Such errors include
issues with the sandwich variance estimation or Jackknife variance estimation
that led to non-positive definite variance matrices. We also considered bias-
eliminated CI coverage to adjust for the impact of bias in the coverage
results morris_using_2019. This involved calculating the frequency at which a
CI generated from each of the 1000 simulated datasets contained the sample
mean of the 1000 MRD estimates, for each simulation scenario.

In Morris et al morris_using_2019 they state 4 potential reasons for CI
undercoverage when analysing simulation results: (i) MRD estimation bias ,
(ii) the MRD distribution sampling methods for constructing the CIs
underestimate the empirical standard deviation of the MRD sampling
distribution, (iii) the distribution of the MRD is not normal and intervals
have been constructed assuming normality (this only concerns the Jackknife CI
using approach 1 in Section 4.4.1, or (iv) the standard error estimate of the
MRD distribution for a given CI construction method is too variable. They also
note that over-coverage can occur when the square of the standard error
estimate of the MRD distribution for a given CI construction method is larger
than the square of the empirical standard deviation of the MRD sampling
distribution. Therefore, for each simulation scenario, we calculated the
empirical bias of the MRD estimates at each trial visit to check point (i) of
CI undercoverage. We collected the standard deviation of the MRD sample from
each CI construction method as standard error estimates $\widehat{SE}$ of the
MRD sampling distribution from that CI method for each simulation, and
calculated, for each CI construction method,
$\frac{1000^{-1}\sum\widehat{SE}}{{SD}_{\widehat{MRD}}}$, the ratio of mean
standard error estimate to the standard deviation of the MRD empirical
distribution ${SD}_{\widehat{MRD}}$, where the latter was calculated by taking
the standard deviation of the 1000 MRD point estimates obtained from the
simulations. This ratio would allow us to assess point (ii) and (iv) of
potential CI undercoverage or potential CI overcoverage.

Additionally, we recorded the number of times that the compared methods
encountered any errors in CI construction, that is, construction failures.
Following Morris et al. (2019) morris_using_2019, we also calculated the Monte
Carlo standard errors for the empirical coverage as a measure of simulation
uncertainty due to using a finite number of simulations (in our case, $1000$
Monte Carlo simulations):

$\displaystyle\mbox{SE}_{\mbox{MC}}=\sqrt{\frac{\widehat{\text{Coverage}}(1-\widehat{\text{Coverage}})}{n_{\mbox{sim}}}},$
(16)

where $\widehat{\text{Coverage}}$ is the empirical coverage and
${n_{\mbox{sim}}}$ is the number of simulations performed
(${n_{\mbox{sim}}}=1000$ in our case).

Finally, we also reported the relative computation time that it took to
construct nonparametric bootstrap-, LEF bootstrap- and Jackknife-based CIs
compared to the sandwich variance estimator-based CIs.

####  5.1.6 Computational resources

The simulations were conducted using the R statistical programming language
(version 4.1.3) r_core_team_r_2023. All R code files are available at
https://github.com/juliettelimozin/Multiple-trial-emulation-IPTW-MSM-CIs.
Point estimation and the CI construction based on the sandwich variance
estimator were performed using the R package TrialEmulation
rezvani_trialemulation_2023 (version 0.0.0.3) and some of the key functions in
the R package are summarised in Section D of the Supplementary Materials. We
used packages doParallel and doRNG to parallelise the bootstrap construction
method over multiple cores in the University of Cambridge high-performance
computing (HPC) cluster. The 1000 simulations for each simulation scenario
took about 5 to 12 hours to complete using 67 cores in the HPC cluster.

###  5.2 Results

The results discussion are organised as follows: we first present the CI
coverage results, then discuss the estimation bias and ratio of mean estimated
standard error to empirical standard deviation of MRD sampling distribution to
explain some of the under or over-coverage results we observed in the
simulation study. We conclude with a discussion on computational time for the
different CI construction methods. Our discussions on Monte Carlo standard
errors and CI construction failure rates can be found in Section E.1 and E.2
of the Supplementary Materials. Since we observed very little difference in
the various performance measures between the two LEF bootstrap methods we will
colloquially refer to them as one when discussing the simulation results. We
refer to the CIs constructed using Approach 1 from Jackknife resampling as
’Jackknife Wald CIs’ and using Approach 2 as ’Jackknife Multivariate Normal
(MVN) CIs’.

####  5.2.1 CI coverage

Figure 1 presents the empirical coverage rates of the CIs based on the
compared methods when the event rate was low, which we will discuss in detail
below. Empirical coverage rates under medium and high event rates can be found
in Figures 13 and 14 of the Supplementary Materials. Figures 19, 20 and 21 of
the Supplementary Materials also presents the bias-eliminated empirical
coverage results of the CIs, which are slightly different from empirical
coverage results for scenarios with non-negligible empirical biases. However,
these differences are not graphically drastic for comparing the CIs methods.

For all methods, we noticed a general pattern of increasing coverage to the
nominal 95% level as the sample size increased, which is not surprising as
more information was available to construct the CIs. There is also a general
trend of all CIs’ coverage decreasing at later visits. This is likely due to
data sparsity at later visits making estimation and inference inaccurate, as
we can see in Table 12 of the Supplementary Materials that we saw very little
to no events at later trial visits in most data scenarios, even when the
sample size was large ($n=5000$).

_Impact of treatment prevalence and sample size._ In scenarios with low/medium
treatment prevalence ($\alpha_{a}=-1,0$) and small/medium sample sizes
($n=200,1000$), nonparametric bootstrap-based CIs consistently exhibited lower
coverage compared to LEF bootstrap-based CIs. For large samples ($n=5000$),
the bootstrap methods performed similarly. This discrepancy may stem from the
random nature of bootstrap resampling, leading to extreme cases with high
treatment imbalance or very few patients with outcome events when sample sizes
were small/medium and event rate was low, thereby affecting the stability of
MSM parameter estimates and CI coverage. Conversely, LEF bootstrapping yielded
more reliable CIs as it avoided the need for matrix inversion and parameter
re-estimation in each bootstrap sample. The non-nominal coverage of
nonparametric bootstrap-based CIs for small/medium sample sizes is also
consistent with the findings of Austin (2022) austin_bootstrap_2022 in the
context of comparing bootstrap standard error to asymptotic variance
estimators when using IPTWs among other weighting methods
austin_bootstrap_2022.

Furthermore, the sandwich variance estimator-based CIs showed poor coverage in
the small sample size scenarios ($n=200$) with low or medium treatment
prevalence ($\alpha_{a}=-1,0$), never achieving the nominal 95% level and
sometimes falling as low as 70%. This poor performance might be explained by
1) the sandwich variance estimation suffering from small sample bias 2) the
large number of construction failures for the sandwich variance estimator-
based CIs when the sample sizes were small ($n=200$); see Section E.2 of the
Supplementary Materials for more details. These results also align with
previous findings on the limitations of the sandwich variance estimator in
small sample cases. Some methods have been proposed to amend this issue
pan_small-sample_2002; rogers_assessment_2021.

For low treatment prevalence ($\alpha_{a}=-1$) and small sample size ($n=200$)
scenarios, the Jackknife Wald CIs consistently maintained the nominal 95%
coverage at all visits, however this coverage decreased as the treatment
prevalence increased ($\alpha_{a}=0$ or $1$). Similar trends could be observed
for the Jackknife MVN CIs, although at low treatment prevalence their coverage
was slightly under the nominal rate.

High treatment prevalence ($\alpha_{a}=1$) resulted in considerable
disparities in CI coverage across the methods. Coverage decayed considerably
at later visits across most methods, with nominal coverage rarely achieved at
post-baseline visit ($k\geq 1$). Nonparametric bootstrap-based CIs had the
most prominent decrease of coverage over time, while the sandwich variance
estimator-based CIs, though not reaching nominal coverage, performed better
than nonparametric bootstrap-based CIs. For small sample sizes, both
Jackknife-based CIs exhibited very poor coverage that slightly increased at
later visits. With small and medium sample sizes, the LEF bootstrap-based CIs
performed similarly as the sandwich variance estimator-based CIs and much
better than nonparametric bootstrap-based CIs. With large sample sizes,
sandwich variance estimator-based CIs approached nominal coverage at earlier
visits but remained suboptimal (around 70% coverage) at later visits, while
all bootstrap methods lagged by about 10%. This is the only simulation
scenario in Figure 1 where the coverage of the sandwich variance estimator-
based CIs surpassed that of LEF bootstrap-based CIs. This reversal in
performance highlights the impact of sample size on the performance of CI
construction methods. With small sample sizes and low event rates, the
sandwich variance estimator exhibited finite-sample bias and encountered many
construction failures. This resulted in poorer CI coverage than the LEF
bootstrap methods. With larger sample sizes, the issues of finite-sample bias
and construction failures for the sandwich variance estimator diminished. On
the other hand, high treatment prevalence and low event rates exacerbated the
treatment imbalances at later visits as seen in the data tabulation examples
in Tables 10, 11 and 12 in the Supplementary Materials, which greatly affected
the performance of the LEF bootstrap and Jackknife methods, even if the sample
sizes were large.

However, caution must be taken when examining the coverage results for the
sandwich variance estimator-based CIs. Indeed Austin (2016)
austin_variance_2016 has shown that this estimator provides a conservative
estimate of the variance matrix, leading to wider and inefficient CIs. This
might have compensated the undercoverage at later visits due to larger
estimation variability and therefore the sandwich variance estimator-based CIs
may be more likely to contain the true MRD.

_Impact of confounding strength and outcome event rates._ We observed minimal
impact on CIs’ coverage with varying levels of confounding strength.

Increasing outcome event rates improved the CI coverage for all methods, as
shown in Figures 13 and 14 in the Supplementary Materials. The impact of
treatment prevalence on Jackknife Wald CIs coverage was reduced when
increasing the outcome event rates, as Jackknife Wald CIs maintained nominal
coverage rates at low/medium treatment prevalence scenarios
($\alpha_{a}=-1,0$) with small sample sizes ($n=200$). Jackknife MVN CIs were
exhibiting over-coverage in such scenarios when the outcome event rate was
high.

With medium event rates and small sample sizes, the sandwich variance
estimator-based CIs and non-parametric bootstrap-based CIs had better coverage
than when the event rate was low, but still performed worse than the LEF
bootstrap-based CIs. This is expected because the sandwich variance estimator
was likely less biased and its construction failures were much less frequent
under higher event rates. Likewise, it was less probable to draw bootstrap
samples with very few or no events, which helped to improve the performance of
the non-parametric bootstrap-based CIs. The patterns of the CIs methods were
similar to those under low event rates when we examined the scenarios with
medium and large sample sizes.

With high event rates and small sample sizes, the performance of the sandwich
variance estimator-based CIs was further improved, but not always better than
the LEF bootstrap-based CIs. But compared with non-parametric bootstrap-based
CIs, the LEF bootstrap-based CIs still performed better under small sample
sizes. With high event rates and medium/large sample sizes, all methods
performed similarly.

Figure 1: Coverage of the CIs under low event rates. Bootstrap: CIs
constructed by nonparametric bootstrap; LEF both: CIs constructed by applying
Approach 2 of LEF bootstrap; LEF outcome: CIs constructed by applying Approach
1 of LEF bootstrap; Jackknife Wald: CIs constructed by applying Approach 1 of
Jackknife resampling; Jackknife MVN: CIs constructed by applying Approach 2 of
Jackknife resampling; Sandwich: CIs based on the sandwich variance estimator

####  5.2.2 Empirical bias of the MRD estimates

Figure 2 presented the empirical biases of the MRD estimates under various
scenarios. We found minimal biases at earlier visits ($k<2$) across various
scenarios of sample size, confounding strength, treatment prevalence and event
rate. But with small sample sizes ($n=200$), the MRD estimates at later visits
($k>3$) exhibited non-negligible biases, especially when the treatment
prevalence was low ($\alpha_{a}=-1$) or high ($\alpha_{a}=1$) and confounding
was the strongest ($\alpha_{c}=0.9$). These are likely finite sample biases,
since, as expected, they approached zero consistently as sample sizes
increased. The increasing biases at later visits coincide with the
deterioration in CI coverage that we observed across all data scenarios at
later visits (point (i) of sources of undercoverage according to Morris et al.
morris_using_2019). Table 12 of the Supplementary Materials provides an
example of the number of events stratified by treatment arm and trial visit
that we observed from one simulated data set in each simulation scenario when
the event rate was low, and we can observe that very few events happened at
later trial visits in these examples, even when the sample size was large
($n=5000$).

The biases under small sample sizes showed different directions across the
scenarios of treatment prevalence and confounding strength. For example, in
scenarios with high treatment prevalence, the biases were negative, which
indicated treatment effect overestimation at later visits. This is likely
because, when treatment reduced the risk of the outcome event and patients
were more likely to receive the treatment than remain untreated due to the
high prevalence of treatment, there were more at-risk patients in the treated
group than in the untreated group at later visits. As a result, due to reduced
sample sizes, the event incidence in the untreated group can be overestimated
at later visits, which led to the overestimation of the treatment effect. We
can observe this phenomenon in Table 10 where with high treatment prevalence,
we observed more patients in the treatment arm but the proportion of events in
the treatment arm was smaller than the proportion of events in the control
arm. Conversely, scenarios with low treatment prevalence showed positive
biases at later visits, indicating treatment effect underestimation. The
positive biases seem to be exacerbated by stronger confounding, while the
negative biases were mitigated. This is not surprising because when strong
confounding was not adequately adjusted for by IPTW under small sample sizes,
treated patients would be at higher risks of the event, which would increase
the positive biases but reduce the negative biases at later visits.

The patterns of the biases under medium event rate scenarios resembled those
under low event rates but with smaller negative biases under high treatment
prevalence and larger positive biases under low treatment prevalence. Under
high event rate scenarios, strong confounding further shifted biases from
negative to positive. This suggests that, when sample sizes were small, the
impact of strong confounding on bias outweighed that of treatment prevalence
when event rates were higher.

We note that patterns of empirical standard deviations of MRD estimation
largely followed those of empirical biases across treatment prevalence and
confounding strength scenarios and can be found in Figure 8. That is, the
scenarios with larger biases also exhibited larger standard deviations. As
expected, larger sample sizes mitigated this variability by providing more
information. Results for the root-Mean Squared Error (MSE) of MRD estimation
were consistent with the results for bias and empirical standard deviations
and can be found in Figure 9 in Supplementary Materials.

Figure 2: Empirical biases of the MRD estimates in various simulation
scenarios.

####  5.2.3 Ratio of mean estimated standard error to empirical standard
deviation of sampling distribution

Figures 3 and 4 present the ratio of mean estimated standard error to
empirical standard deviation of sampling distribution of the MRD for scenarios
with low event rates, with the minimum and maximum values of the bar around
the means representing the 1st and 3rd quantile of this ratio respectively, to
display the variability of the standard error estimation. The ratio figures
for medium and high event rate scenarios can be found in Figures 15, 16, 18
and 18 of the Supplementary Materials.

In general, with nonparametric bootstrap-based CI construction we observed a
ratio that was on average consistently below 1 at all visits, meaning
nonparametric bootstrap resampling tended to underestimate the standard
deviation of the sampling distribution of the MRD. The standard deviation
estimation also displayed more variability at later visits with the quantile
bars of the ratios becoming longer. This is consistent with the coverage
results from the previous section, where we would expect the nonparametric
bootstrap-based CIs to have poor coverage if the sampling variability of the
MRD sampling distribution was underestimated and the estimates were highly
variable morris_using_2019. In contrast, for small sample size scenarios
($n=200$) the Jackknife standard error estimate of the MRD that was used for
constructing the Jackknife Wald CIs was more accurately estimating the
empirical standard deviation of the MRD sampling distribution at low or medium
treatment prevalence ($\alpha_{a}=-1,0$) as the ratio was close to 1
consistently, with small variability in this standard deviation estimation,
which coincides with where we observed nominal or almost nominal coverage
rates of these CIs. However the ratio for Jackknife MVN CIs is much larger
than 1 for all scenarios with small sample sizes, which may explain the
overcoverage results we observed in the CI coverage results for some
scenarios.

Both LEF bootstrap sampling of the MRD and the multivariate normal sampling of
the MRD using the sandwich variance estimator resulted in both over- and
underestimation of the empirical standard deviation of the MRD. For small
sample sizes ($n=200$), the ratios for both methods was largely above 1 then
gradually fell below 1. When the treatment prevalence was high
($\alpha_{a}=1$), their ratio was much smaller than that of nonparametric
bootstrap at the later visits. This larger underestimation of the sampling
variability for later visits compared to nonparametric bootstrap may be
counter-intuitive given that LEF bootstrap- and sandwich variance estimator-
based CIs had higher coverage than nonparametric bootstrap-based CIs, however
the overestimation of standard deviation for LEF bootstrap- and sandwich
variance estimator-based CI construction may have some protection for better
coverage at earlier visit, which may also compensate the underestimation at
later visit in term of CI coverage. The LEF bootstrap- and sandwich variance
estimator-based standard error ratios were also much less variable than the
nonparametric bootstrap sampling’s ratios which may explain why LEF bootstrap-
and sandwich variance estimator-based CI undercoverage was not as poor as
nonparametric bootstrap-based CIs.

In medium and large sample sizes ($n=1000,5000$) the LEF bootstrap sampling
and the multivariate normal sampling of the MRD using the sandwich variance
estimator tended to overestimate the empirical standard deviation of the MRD,
with the ratio being above 1 at most visits, which could explain some of the
CI overcoverage we observed. The increase of sample size also reduces under or
overestimation of the empirical standard deviation across all methods, as the
ratios approach 1 as sample size increases, and the same can be observed when
increasing outcome event rates. This is consistent with the coverage results
we observed earlier, where an increase of sample size and outcome event rate
resulted in the various CI construction methods converging to similar, closer
to nominal coverage. However in large sample sizes and low or high treatment
prevalence, we can still observe ratios to be lower than one at later visits
and that are also highly variable. The underestimation and high variability of
the standard errors at later visits could explain the lower coverage we
observed at later visits for all CI methods in the previous section, even at
large sample sizes.

Similarly to the empirical coverage results, the confounding strength does not
seem to have affected sampling variability estimation.

Figure 3: Ratio of mean estimated standard error to empirical standard
deviation of MRD sampling distribution under low event rates. Bootstrap: CIs
constructed by nonparametric bootstrap; LEF both: CIs constructed by applying
Approach 2 of LEF bootstrap; LEF outcome: CIs constructed by applying Approach
1 of LEF bootstrap; Jackknife Wald: CIs constructed by applying Approach 1 of
Jackknife resampling; Sandwich: CIs based on the sandwich variance estimator

Figure 4: Ratio of mean estimated standard error to empirical standard
deviation of MRD sampling distribution under low event rates. Jackknife MVN:
CIs constructed by applying Approach 2 of Jackknife resampling

####  5.2.4 Computation time

Table 4 summarises the average computation time for constructing a CI of the
MRD using each method, relative to the time for constructing the sandwich
variance estimator-based CI. Each average is over the 1000 datasets per
scenario and over the various scenarios of treatment prevalence and
confounding strength. On average, nonparametric bootstrap CIs took about 3 to
8 times longer to compute than the sandwich variance estimator-based CIs. The
computation time increases at a fast rate as sample sizes increase. LEF
bootstrap-based CIs, using both Approaches 1 and 2, took about 1.6 to 2.5
times longer to compute compared to the sandwich variance estimator-based CIs.
LEF bootstrap CIs are not as affected by sample size in terms of computation
time compared to nonparametric bootstrap-based CIs.

Jackknife Wald CIs had very similar computation times to LEF bootstrap-based
CIs at small sample sizes ($n=200$) whereas Jackknife MVN CIs took 3 times
longer to construct compared to sandwich variance estimator-based CIs, most
likely because this method requires two sampling tasks: the Jackknife
resampling process and the MVN sampling process.

The large gain in computational efficiency from nonparametric to LEF bootstrap
is very much of interest; however, the sandwich variance estimator-based
method remains the fastest for CI construction.

Table 4: Summary of average computation time of the methods for constructing CIs in simulations with the sandwich variance estimator-based CIs as the reference. Bootstrap: CIs constructed by nonparametric bootstrap; LEF both: CIs constructed by applying Approach 2 of LEF bootstrap; LEF outcome: CIs constructed by applying Approach 1 of LEF bootstrap; Jackknife Wald: CIs constructed by applying Approach 1 of Jackknife resampling; Jackknife MVN: CIs constructed by applying Approach 2 of Jackknife resampling; Sandwich: CIs based on the sandwich variance estimator  Outcome even rate | Sample size | Bootstrap | LEF | LEF | Jackknife | Jackknife | Sandwich
---|---|---|---|---|---|---|---
event rate |  |  | outcome | both | Wald | MVN |
Low | 200 | 2.94 | 1.83 | 1.62 | 1.78 | 3.00 | 1
| 1000 | 4.28 | 1.92 | 1.74 |  |  | 1
| 5000 | 8.39 | 2.66 | 2.50 |  |  | 1
Medium | 200 | 2.99 | 1.88 | 1.66 | 1.81 | 3.04 | 1
| 1000 | 4.31 | 1.99 | 1.77 |  |  | 1
| 5000 | 8.39 | 2.68 | 2.51 |  |  | 1
High | 200 | 2.97 | 1.88 | 1.66 | 1.80 | 3.02 | 1
| 1000 | 4.23 | 1.96 | 1.77 |  |  | 1
| 5000 | 7.99 | 2.65 | 2.48 |  |  | 1

####  5.2.5 Summary

The results of our simulation study suggest that, in comparison to the
sandwich variance estimator-based CIs, LEF bootstrap-based CIs provided better
coverage in settings with small/moderate sample sizes. The performance of the
nonparametric bootstrap-based CIs was largely affected by large treatment
imbalances, especially in settings with small/moderate sample sizes and low
event rates, while the LEF bootstrap-based method had better performance and
was computationally more efficient in such settings. For small sample sizes
and low treatment prevalence data scenarios, Jackknife Wald CIs achieved
nominal coverage but were also sensitive to treatment imbalances. All methods
performed poorly when the prevalence of treatment was high. Since the use of
sequential trials is particularly useful for data scenarios with small numbers
of patients initiating treatments at any given time (low treatment prevalence)
and with low event rates , LEF bootstrap could be an advantageous alternative
to the sandwich variance estimator and the nonparametric bootstrap in
sequential emulated trials with small/moderate sample sizes. Jackknife Wald
CIs also performed well in terms of coverage in such scenarios however if
computational efficiency is an issue we would not be recommended this
construction method when the number of patient exceeds the number of bootstrap
samples one would generate for the LEF bootstrap as the Jackknife resampling
would become computationally inefficient compared to LEF bootstrap. Jackknife
MVN CIs encountered too many construction failures during the simulations as
seen in Section E.2 so it would be unadvisable to use this method. With large
sample sizes, we recommend sandwich variance estimator-based CIs as they
performed better than two bootstrap-based CIs and were also computationally
efficient.

The results of our simulation study largely highlighted how different data
scenarios affect the performance of CIs and how one needs to carefully
consider the longitudinal data’s structure (sample size, outcome even rate,
treatment prevalence) to make a decision on which CI construction method one
should choose.

##  6 Application to the HERS data

In this section, we applied the methods described in previous sections to the
HERS data. Following the protocol in Table 1 and using the 584 women assessed
at visit 8 of the HERS, we emulated 5 trials $m=0,ldots,4$ with the trials
starting at visits 8, 9, 10, 11 and 12, respectively. Patients should have no
history of receiving HAART before the baseline of each trial in order to be
eligible for enrollment in that trial. Dependent censoring was present with
175 women lost to follow-up from visit 8 to visit 12. Of the 584 women
observed at visit 8, 17 women were lost to follow-up after visit 8, 24 women
after visit 9, 24 women after visit 10 and 110 women after visit 11. 385 women
were observed up to visit 12. A total of 24 events were recorded, with 4 women
experiencing the event between visit 8 and 9, 7 women between visit 9 and 10,
one woman between visit 10 and 11, 2 women between visit 11 and 12 and 10
women after visit 12. Additionally, 491 women were eligible to trial 0, 401
women to trial 1, 320 women to trial 2, 255 to trial 3 and 179 trial 4,
resulting in 1646 patient observations by pooling the 5 emulated trials. A
summary of the data tabulation by treatment group, outcome event and dependent
censoring indicator at each emulated trial is given in Table 5.

Table 5: Data tabulation of the HERS data after expanding the dataset into sequentially emulated trials. The numbers represent the aggregated observations over the trial visits. Treatment: treatment protocol indicator, 0: never treated, 1: always treated. Outcome: event indicator. C: dependent censoring due to loss to follow-up indicator. |  |  Trial index $m$ | 0 | 1 | 2 | 3 | 4
---|---|---|---|---|---|---|---
Treatment | Outcome | $C$ |  |  |  |  |
0 | 0 | 0 | 390 | 315 | 250 | 175 | 155
|  | 1 | 14 | 16 | 15 | 49 | 0
| 1 | 0 | 11 | 8 | 5 | 5 | 4
1 | 0 | 0 | 73 | 54 | 49 | 25 | 19
|  | 1 | 1 | 5 | 0 | 1 | 0
| 1 | 0 | 2 | 3 | 1 | 0 | 1
|  | Total eligible in trial | 491 | 401 | 320 | 255 | 179

The treatment process for the denominator term in formula (5) was modelled by
two logistic regressions stratified by treatment received at the immediately
previous visit $A_{t_{j-1}}$ ($j=1,\ldots,4$). As HAART treatment was a self-
reported status between the last visit and the current visit, we included the
following covariates: CD4 count (after square root transformation and
standardization) and HIV viral load (after $\log_{10}$ transformation and
standardization) measured at the previous two visits, HIV symptoms at the
previous visit, ethnicity (Caucasian, Black, other) and study sites. The
numerator term in formula (5) was estimated by the marginal probabilities of
receiving HAART stratified by previous treatments $A_{t_{j-1}}$.

Similarly, the censoring process for the denominator term of the formula (6)
was modelled by two logistic regressions stratified by the previous treatment
$A_{t_{j-1}}$. The following covariates were included: CD4 count, HIV viral
load and HIV symptoms measured at the current and previous visits, ethnicity
and study sites. The numerator term in formula (6) was estimated by the
marginal probabilities of remaining in the HERS cohort, stratified by the
previous treatment $A_{t_{j-1}}$.

We specified the MSM for the emulated trials in terms of the discrete-time
hazard of all-cause mortality by applying a pooled logistic regression with
the following covariate terms included: the intercept term, the main treatment
effect, the main effects of CD4 count and HIV viral load at previous two
visits before the trial baseline, the main effects of ethnicity and study
sites, and the interaction between treatment and CD4 count at the visit before
trial baseline. A summary of the fitted MSM is provided in Table 6.

Table 6: Results for the fitted MSM for the sequentially emulated trials using the HERS data. Robust standard error: standard error estimate from the robust sandwich variance estimate of the MSM. | Estimate | Robust standard error
---|---|---
Intercept | -5.159 | 1.131
Assigned treatment: HAART vs. non-HAART | -0.015 | 0.549
CD4 count at 1 visit before trial baseline | 0.146 | 0.380
CD4 count at 2 visits before trial baseline | -0.277 | 0.372
Viral load at 1 visit before trial baseline | 0.015 | 0.354
Viral load at 2 visits before trial baseline | 0.615 | 0.275
Site 2 vs Site 1 | 0.503 | 1.148
Site 3 vs Site 1 | 0.071 | 0.936
Caucasian vs. Black | -0.182 | 0.757
Other ethnicity vs. Black | 0.977 | 1.078
Interaction between assigned treatment and | -0.004 | 0.554
CD4 count at 1 visit before trial baseline |  |

CIs of the MRD were constructed as described in Section 4 with 500 bootstrap
samples for nonparametric and LEF bootstrap, and drawing a sample of size
$S=500$ from the multivariate normal distribution in the sandwich variance
estimator method for constructing CIs.

Figure 5 presents the estimated MRD, which is the difference in risk of all-
cause mortality when all patients eligible to trial 0 adhere to the treatment
strategies assigned ( HAART vs. no HAART). Figure 5 also includes different
95% CIs of the MRD obtained by the four methods described in Section 4. We
note that the results are not statistically significant given that all 4 CIs
cover zero.

Figure 5:  Estimates and 95% CIs of the MRD of all-cause mortality under HAART
treatment strategy and the control strategy for patients enrolled in the first
emulated trial of the HERS data. Nonparametric Boot.: CIs constructed by
nonparametric bootstrap; LEF both: CIs constructed by applying Approach 2 of
LEF bootstrap; LEF outcome: CIs constructed by applying Approach 1 of LEF
bootstrap; Sandwich: CIs based on the sandwich variance estimator

In addition, we wanted to compare the per-protocol effect estimation and
inference using a more common approach in practice clare_causal_2019, that is,
fitting a standard MSM using IPW. This approach was first introduced by Robins
et al. robins_marginal_2000 and has been applied to survival outcomes by
Hernán et al. hernan_marginal_2000. This approach, which we will denote by
MSM-IPW keogh_causal_2023, does not require artificial censoring of patients’
follow-up at the point of treatment non-adherence, and instead, it consists of
estimating the inverse probability of treatment weights (and additionally,
inverse probability of censoring weights if dependent censoring is present)
and the MSM for the potential outcomes using the full observational cohort
data that have not been artificially censored. To apply the MSM-IPW approach
to the subcohort of 584 HERS patients, the inverse probability of treatment
and censoring weights to adjust for time-varying confounding of the treatment
sequence and dependent censoring were estimated via logistic regressions, with
the same model specification as for the analysis of the sequentially emulated
trial described earlier. A weighted pooled logistic regression was fitted to
the full observational data in the HERS subcohort to estimate the MSM
parameters. A summary of this fitted MSM is provided in Table 7. Similarly to
the analysis of the sequential trials, the model specification included an
intercept term, the main effect of cumulative exposure to HAART, the main
effect of CD4 count and viral load at 2 previous visits before cohort baseline
(visit 8), main effects of ethnicity and study sites and interaction between
cumulative exposure to HAART and CD4 count at the previous visit before cohort
baseline.

Table 7: Results for the fitted MSM in the MSM-IPW approach using the HERS data. Std. Error: model-based standard error estimate of the model coefficients. | Estimate | Std. Error | z value | Pr($>$$|$z$|$)
---|---|---|---|---
Intercept | -5.787 | 0.790 | -7.329 | 0.000
Cumulative exposure to HAART | 0.001 | 0.002 | 0.409 | 0.682
CD4 count at 1 visit before cohort baseline | -0.141 | 0.666 | -0.211 | 0.833
CD4 count at 2 visits before cohort baseline | 0.033 | 0.519 | 0.064 | 0.949
HIV viral load at 1 visit before cohort baseline | -0.252 | 0.333 | -0.758 | 0.449
HIV viral load at 1 visit before cohort baseline | 0.838 | 0.368 | 2.276 | 0.023
Site 2 vs Site 1 | 0.669 | 1.039 | 0.643 | 0.520
Site 3 vs Site 1 | -0.109 | 1.462 | -0.075 | 0.941
Caucasian vs Black | -0.151 | 0.646 | -0.234 | 0.815
Other ethnicity vs Black | 0.487 | 0.700 | 0.696 | 0.486
Interaction between cumulative exposure to HAART and | -0.001 | 0.001 | -0.674 | 0.500
CD4 count at 1 visit before cohort baseline |  |  |  |

The MRD was then estimated by predicting the discrete-time hazards from this
fitted MSM on two datasets containing all patients in the HERS subcohort with
their cohort baseline covariate information, with one dataset as having them
all be ‘always-treated’ and the other dataset having them as all ‘never-
treated’, in a very similar manner to how the MRD is estimated from the MSM
fitted in the sequential trial emulation approach We then constructed 95%
pivot bootstrap CIs of the MRDs using 500 nonparametric bootstrap samples,
similarly to the sequential trial emulation approach.

Figure 6 presents the MRD estimates from the sequential trial emulation
approach and the MSM-IPW approach. The CIs presented in Figure 6 were both
generated via nonparametric bootstrap. We can note that although the MRD
estimates are very similar in both approaches, the bootstrap CIs obtained from
the MSM-IPW approach are far wider than those from the sequential trial
emulation approach. This suggests that estimation via MSM-IPW might be less
efficient than estimation via sequential trial emulation, which might be
attributed to the larger variability of inverse probability weights in the
traditional MSM-IPW analysis. Table 8 provides some summary statistics of the
estimated weights in both approaches. It is clear that the weights from the
MSM-IPW approach were more extreme and more variable than the weights in the
sequential trial emulation approach. These results are in line with the
simulation findings of Keogh et al. keogh_causal_2023 who compared the
performance of treatment effect estimation under both approaches when the MSMs
for counterfactual hazards were based on additive hazards models and found
that the relative efficiency of the sequential trial emulation approach is
related to its less extreme weights compared with those weights from the MSM-
IPW approach.

Table 8: Summary statistics of the estimated weights from the sequential trial emulation and traditional MSM-IPW approaches in the HERS data analysis. Statistic | Weights in sequential trial emulation | Weights in the MSM-IPW approach
---|---|---
Mean | 0.99 | 1.15
Standard deviation | 0.16 | 0.76
Minimum | 0.50 | 0.20
Maximum | 2.53 | 23.73
1st percentile | 0.65 | 0.48
99th percentile | 1.65 | 3.24
Figure 6: Marginal risk difference estimates and $95\%$ CIs based on
nonparametric bootstrap when using the sequential trial emulation and
traditional MSM-IPW approaches.

##  7 Discussion

In this paper, we focused on the application of sequential trial emulation to
estimate and make inferences about the per-protocol effect of treatments on a
survival outcome in terms of counterfactual MRDs over time. Motivated by its
computational efficiency, we extended the LEF bootstrap to construct CIs for
MRDs and compared its performance with four other CI methods based on the
nonparametric bootstrap, the sandwich variance estimator and Jackknife
resampling using simulations. Based on our simulation results, we recommend
sandwich variance estimator-based CIs for sequential trial emulation with
large sample sizes and LEF bootstrap-based CIs for settings with
small/moderate sample sizes.

Although dependent censoring was not included in the data generating mechanism
of our simulation study, we expect that the relative performance of CIs
constructed using LEF bootstrap and the sandwich variance estimator that we
found in our simulation study would also hold in settings where inverse
probability of censoring weighting is additionally used to handle dependent
censoring. Inference would still face the same issues, including additional
uncertainty caused by also having to estimate the inverse probability of
censoring weights and the finite-sample bias that may be exacerbated by loss
to follow-up.

We primarily focused on pooled logistic regressions for fitting MSMs and did
not explore alternative survival time models, such as additive hazards models
keogh_causal_2023. Depending on the scenarios in which the models are
suitable, different outcome modelling strategies could have been considered.
For instance, additive hazard models can be used in moderate to high event
rates settings, while in low event settings, they might result in negative
hazard estimates.

We also note that the analysis in this paper focused on parametric models for
IPW weight estimation. Recent developments in nonparametric data-adaptive
models have made them widely used in biomedical research. This is partly
because the implementation of these models in programming languages has become
more automated, making them more accessible. Moreover, data-adaptive models
are appealing because they provide flexibility when specifying a model due to
their nonparametric nature. Because the consistency of inverse probability
weighted estimators of treatment effects relies on the correct specification
of treatment and censoring models, data-adaptive models may be the more
judicious choice for model specification when we are performing target trial
emulation with complex relationships between treatment, confounders, outcome
and censoring processes. However, the subsequent inference can be challenging
schuler_targeted_2017; chernozhukov_doubledebiased_2016;
cai_nonparametric_2020.

Additionally, throughout this paper, it was assumed that there were no
unmeasured confounders. However, in practice, this assumption is unlikely to
hold as it is challenging to identify and measure all potential confounders in
observational databases. For unmeasured confounding, instrumental variables
and sensitivity analysis are commonly employed approaches. However, applying
these methods in longitudinal settings poses additional difficulties.
Recently, Tan (2023) extended the marginal sensitivity models to longitudinal
settings to perform sensitivity analysis for the sequential unmeasured
confounding assumption tan_sensitivity_2023, which might shed light on the
development of sensitivity analysis approaches for unverifiable assumption in
target trial emulation.

##  Appendix A Sequential Trial Emulation illustration

S

Figure 7: An illustration of how a patient’s observed data is utilised in
sequentially emulated trials.

Three trials, starting at visits at time $t=0,1,2$, are sequentially emulated.
Patient 1 never initiated treatment for the duration of the study until they
were lost to follow-up after the second visit, so they are eligible as a
control in all three trials. Patient 2 initiated treatment (marked by the pink
highlight) at time 0 and kept taking the treatment until the end of follow-up
so they are in the treated group for the first trial, but not eligible for the
second and third trial under the assumption that patients are only eligible if
there is no prior history of treatment. Patient 3 initiated treatment at time
0 but stopped from time 1, so they are artificially censored from $t=1$. They
are included in the treatment group for the first trial, but their observed
data is not included in other trials because they had a prior history of
treatment. Patient 4 was not taking treatment at time 0 but initiated
treatment from time $t=1$ onwards: hence they are eligible as a control in the
first trial and artificially censored at the time they started treatment, and
they are eligible as a treated in the second trial.

##  Appendix B Causal Assumptions

No Interference.

Patient $i$’s received treatment has no effect on patient $j$’s counterfactual
outcomes ($i\neq j$) daniel_methods_2013, i.e., a patient’s event time is not
affected by other patient’s treatments.

Consistency.

$\forall m,\>k,\>\overline{Y}_{m,k}=\overline{Y}_{m,k}^{\overline{A}_{m,k}}$
and $\overline{L}_{m,k}=\overline{L}_{m,k}^{\overline{A}_{m,k}}$, meaning that
the observed outcomes $\overline{Y}_{m,k}$ and covariates $\overline{L}_{m,k}$
are equal to their counterfactual outcomes and covariates, under the treatment
assignment which they actually received, for every trial visit $k$ in trial
$m$ daniel_methods_2013; young_simulation_2014.

Positivity.

$\forall
m,\>k,\>\Pr[A_{m,k}=a|\overline{A}_{m,k-1}=\overline{a},V,\overline{L}_{m,k},Y_{m,k-1}=0,C_{m,k-1}=0]>0$
and
$\Pr[C_{m,k}=1|C_{m,k-1}=0,\overline{A}_{m,k}=\overline{a},V,\overline{L}_{m,k},Y_{m,k-1}=0]>0$,
i.e., a patient could have been assigned to either treatment or lost to
follow-up at all time conditional on treatment history $a=0,1$ and any history
of $L$ young_simulation_2014.

Sequentially ignorable treatment assignment.

$\forall
m,\>k,\>Y_{m,k}^{\overline{a}},\dots,Y_{m,K}^{\overline{a}}{\,\perp\\!\\!\\!\perp\,}A_{m,k}|V,\overline{L}_{m,k},\overline{A}_{m,k-1},Y_{m,k-1}=0,C_{m,k-1}=0$
for $a=0,1$, i.e., at a given time, conditional on past treatment assignment
and covariate history, there is no unmeasured confounding between future risk,
the current treatment assignment and the current covariates
daniel_methods_2013.

Sequentially ignorable loss to follow up.

$\forall
m,\>k,\>\overline{C}_{m,k}{\,\perp\\!\\!\\!\perp\,}Y_{m,k+1}^{\overline{a}_{k}},\dots,Y_{m,K}^{\overline{a}_{k}}|\overline{C}_{m,k-1}=0$,
$Y_{m,k}=0$, $\overline{A}_{m,k}=\overline{a}$, $V$, $\overline{L}_{m,k}$, for
$a=0,1$ daniel_methods_2013. In other words, at a given time, a patient’s
probability of being lost to follow-up does not depend on their future risk of
event, conditional on the treatment strategy and covariate history up to that
time.

##  Appendix C Pseudo code for simulation study

Algorithm 1 presents a pseudo-code of the Monte Carlo simulation studies. The
R code can be found at https://github.com/juliettelimozin/Multiple-trial-
emulation-IPTW-MSM-CIs.

Algorithm 1 Pseudo-code of Monte Carlo simulation algorithm.

for MC iteration $m=1,2,...,1000$ do

for Data scenario $l=1,2,...,81$ do

Data A $\leftarrow$ Generate data with scenario $l$

Select eligible patients

Fit stabilised IPTWs denominator and numerator models on eligible patients’
data up to the first visit they switched treatment after trial baseline.

Data B $\leftarrow$ the original dataset expanded into sequentially emulated
trials and assign patients to eligible trials

for Each trial do

Artificially censor patients when they no longer adhere to trial baseline
treatment protocol

Estimate each patient’s stabilised IPTW based on covariate history in this
trial. At trial baseline, the weights are always one.

end for

Fit weighted pooled logistic regression on Data B

$v\leftarrow$ MRD estimate for trial 0 patients

$[B_{1},B_{2}]\leftarrow$ lower- and upper-bound of nonparametric bootstrap
pivot CI of MRD

$[LEF^{(1)}_{1},LEF^{(1)}_{2}]\leftarrow$ lower- and upper-bound of LEF
bootstrap pivot CI of MRD using approach 1

$[LEF^{(2)}_{1},LEF^{(2)}_{2}]\leftarrow$ lower- and upper-bound of LEF
bootstrap pivot CI of MRD using approach 2

$[S_{1},S_{2}]\leftarrow$ lower- and upper-bound of sandwich variance
estimator-based CI of MRD

end for

end for

return
$v,[B_{1},B_{2}],[LEF^{(1)}_{1},LEF^{(1)}_{2}],[LEF^{(2)}_{1},LEF^{(2)}_{2}],[S_{1},S_{2}]$

##  Appendix D TrialEmulation package: more details on key package functions

  * •

data_preparation: takes longitudinal data and IP models specifications as
input, and returns an expanded dataset with patients’ data spread across their
eligible sequentially emulated trials. This dataset includes a column
indicating the patient’s inverse probability weights as calculated from the
formulae (6) and (5), according to whether we are conducting an intention-to-
treat or per-protocol analysis. The fit of the models required for these
formulae was done pre-expansion, using the variables indicated in the function
input. For per-protocol analysis, the output expanded data has been
artificially censored whenever patients stop adhering to their baseline
treatment protocol in each emulated trial. This function operates steps 1 and
2 of treatment effect estimation in multiple trials setting presented at the
end of Section 3.

  * •

trial_msm: takes an expanded data with inverse probability weights such as one
returned from data_preparation and the outcome model specification as inputs
and produces a fit.glm type output of the MSM with estimated hazard ratios and
the sandwich variance estimator. This function operates step 3 of treatment
effect estimation in a sequentially trial emulation setting. We can use
predict.glm on the output of this function to estimate cumulative hazards from
a given data, which is what we use to carry out step 5 of the analysis.

  * •

initiators: a summary function combining the two functions above with
additional inputs for computational efficiency factors.

The data_preparation function in the package was used to estimate the
stabilised IPTWs from the simulated datasets, and this function outputs the
sequential trials-expanded and artificially censored data as described in 3,
with each patient’s inverse probability weights at each trial’s follow-up
time. The discrete-time hazard MSM for the outcome was fitted using trial_msm.

The robust sandwich variance estimator is provided in the output of trial_msm
and was used to obtain CIs. The bootstrapping process code, both for
nonparametric and LEF CIs, were not included in the TrialEmulation package so
they were created for the purpose of this paper and can be found in the GitHub
repository.

More details on the TrialEmulation package can be found at
https://cran.r-project.org/web/packages/TrialEmulation/TrialEmulation.pdf

##  Appendix E More simulation results

####  E.0.1 Empirical standard deviation and MSE of the MRD estimates

For all scenarios examined, the empirical standard deviations (Figure 8) were
always higher at later visits than at earlier visits. There are a few possible
reasons for this: 1) due to the non-collapsibility of logistic models, we were
not able to specify an exactly correct MSM for the discrete-time hazard.
Instead, we used a rich model by visits, which might have led to larger
variability in MRD estimates at later visits. 2) Only earlier trials had data
for longer follow-ups, so the information was scarce by design. 3) As
mentioned previously, we rely on IPW to estimate the MRD. The limited amount
of observations at later visits paired with the increased imbalance between
the treatment groups likely resulted in the estimated weights for later trials
visits being biased or highly variable. Consequently, the larger variability
of the weights at later visits would contribute to the unstable estimation of
MRD.

The patterns of empirical standard deviations largely followed those of
empirical biases across treatment prevalence and confounding strength
scenarios. That is, the scenarios with larger biases also exhibited larger
standard deviations. Overall, the impact of treatment prevalence on empirical
standard deviations was more prominent than that of confounding strength,
possibly due to the small range we specified for the coefficients of the time-
varying confounder in the treatment and outcome processes. Interestingly, the
empirical standard deviations increased when event rates were higher, when
fixing the sample size, treatment prevalence and confounding strength. This
phenomenon might be explained by that more patients from the untreated group
were depleted from the at-risk sets due to more event occurrences, which
resulted in larger treatment group imbalances. As expected, larger sample
sizes mitigated this variability by providing more information. Results for
the root-MSE were consistent with the results for bias and empirical standard
deviations and can be found in Figure 9 in Supplementary Materials.

The bias and variability of the MRD estimates could partly be explained by the
variability of the estimated IPTWs across scenarios. Tables 13, 14 and 15 in
the Supplementary Materials show summary statistics of the estimated IPTWs for
one simulated dataset from each scenario. It is clear that confounding
strength had a substantial impact on the how large the IPTWs could be (a
measure of right-skewness of the weights), with stronger confounding leading
to larger variability (standard deviations) of the weights.

Figure 8: Empirical standard deviation of the MRD estimates in various
simulation scenarios.

Figure 9: Root-mean-squared error (MSE) of the MRD estimates in various
simulation scenarios.

###  E.1 Monte Carlo standard error

Finally, to assess whether lower-than-nominal coverage in certain data
scenarios may be due to uncertainty induced by Monte Carlo simulation error,
we calculated and illustrated the standard errors of empirical coverage rates
in Figures 10,11 and 12. We observe minimal Monte Carlo standard errors across
all simulated data scenarios, with the error mostly remaining at 1.5%,
especially for small sample sizes where we did not observe much difference in
Monte Carlo error between the CI methods. According to the Monte Carlo
standard error formula for CI coverage, if the expected coverage is 95% and we
want to observe 1.5% standard error or less across the simulations, then
$n_{sim}=\frac{95(100-95)}{1.5^{2}}=211$ simulations should be enough to
observe nominal coverage morris_using_2019. According to this formula, 1000
simulations were therefore sufficient to theoretically achieve nominal
coverage with at most 1.5% error across the Monte Carlo simulations for small
$n=200$ sample size scenarios, where the MC error remained at 1.5%, so the
bias leading to poor coverage is most likely due to the lack of information
for events in one of the treatment arms at later visits.

We also observed lower MC error at earlier follow-ups in large ($n=5000$)
sample size scenarios. For the Monte Carlo errors, there was a clear trend of
increasing at later visits for $n=5000$ scenarios, meaning that the sparsity
of the data was a problem for later visits even for large sample sizes. As the
event rate increases, the MC errors converge to similar values for all CI
methods as well. It is of interest to also note that at large sample sizes,
sandwich variance-based CIs resulted in much lower MC error across the
simulation than the bootstrap CI methods, which may suggest that some of the
lower coverages for the bootstrap CIs at large sample sizes, compared to the
sandwich CIs, may have been due to their performance being affected more by
random errors than the sandwich variance estimator CIs.

Figure 10: Monte Carlo standard error the CIs under low event rates.
Bootstrap: CIs constructed by nonparametric bootstrap; LEF both: CIs
constructed by applying Approach 2 of LEF bootstrap; LEF outcome: CIs
constructed by applying Approach 1 of LEF bootstrap; Jackknife Wald: CIs
constructed by applying Approach 1 of Jackknife resampling; Jackknife MVN: CIs
constructed by applying Approach 2 of Jackknife resampling; Sandwich: CIs
based on the sandwich variance estimator

Figure 11: Monte Carlo standard error the CIs under medium event rates.
Bootstrap: CIs constructed by nonparametric bootstrap; LEF both: CIs
constructed by applying Approach 2 of LEF bootstrap; LEF outcome: CIs
constructed by applying Approach 1 of LEF bootstrap; Jackknife Wald: CIs
constructed by applying Approach 1 of Jackknife resampling; Jackknife MVN: CIs
constructed by applying Approach 2 of Jackknife resampling; Sandwich: CIs
based on the sandwich variance estimator

Figure 12: Monte Carlo standard error the CIs under high event rates.
Bootstrap: CIs constructed by nonparametric bootstrap; LEF both: CIs
constructed by applying Approach 2 of LEF bootstrap; LEF outcome: CIs
constructed by applying Approach 1 of LEF bootstrap; Jackknife Wald: CIs
constructed by applying Approach 1 of Jackknife resampling; Jackknife MVN: CIs
constructed by applying Approach 2 of Jackknife resampling; Sandwich: CIs
based on the sandwich variance estimator

###  E.2 Convergence issues

We observe from Table 9 that whilst we encountered little to no convergence
issues when constructing nonparametric bootstrap- and LEF bootstrap-based CIs
and Jackknife Wald CIs for all data scenarios, there were substantially more
problems with positive definiteness and symmetry when computing the sandwich
variance matrices. This led to a high rate of failures, ranging from 3 to 81%,
in constructing sandwich variance estimator-based CIs as the multivariate
normal distribution for sampling requires a symmetric, positive definite
variance-covariance matrix. There were problems with these matrix properties
for more than half of the simulations when the sample size was $n=200$, whilst
the failure rate for sample size $n=1000$ was between 3% and 37%. Although we
have not included the disaggregation by confounding, we note that the failure
rate for sandwich variance estimator-based CIs was further exacerbated in
scenarios with $n=200$ and high treatment prevalence, averaging at around 81%.
On occasion, we observed that this was due to R being too restrictive with the
matrix properties. In some cases, the sandwich matrix was almost symmetric
with differences between symmetric values at a negligible magnitude of
$10^{-4}$, so discarding such a difference and forcing a symmetric matrix for
sandwich CIs may have been judicious, and may have increased the coverage rate
results for these intervals slightly. For very large sample sizes $n=5000$,
the failure rate was negligible across all methods.

Jackknife MVN CIs encountered further construction issues, where for low and
high treatment prevalence this method failed for roughly half of the
simulations, making this CI construction method unreliable for such data
scenarios.

Table 9: Frequency of CI construction failure aggregated by outcome event
rate, sample size and treatment prevalence. Bootstrap: CIs constructed by
nonparametric bootstrap; LEF both: CIs constructed by applying Approach 2 of
LEF bootstrap; LEF outcome: CIs constructed by applying Approach 1 of LEF
bootstrap; Jackknife Wald: CIs constructed by applying Approach 1 of Jackknife
resampling; Jackknife MVN: CIs constructed by applying Approach 2 of Jackknife
resampling; Sandwich: CIs based on the sandwich variance estimator

Outcome event rate | Sample size | Treatment prevalence | Bootstrap | LEF outcome | LEF both | Sandwich | Jackknife MVN | Jackknife Wald
---|---|---|---|---|---|---|---|---
Low | 200 | -1 | 0.00 | 0.00 | 0.00 | 0.07 | 0.46 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.07 | 0.46 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.07 | 0.46 | 0.00
|  | 0 | 0.00 | 0.00 | 0.00 | 0.27 | 0.03 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.27 | 0.03 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.27 | 0.03 | 0.00
|  | 1 | 0.00 | 0.00 | 0.00 | 0.81 | 0.51 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.81 | 0.51 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.81 | 0.51 | 0.00
| 1000 | -1 | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  | 0 | 0.00 | 0.00 | 0.00 | 0.03 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.03 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.03 |  |
|  | 1 | 0.00 | 0.00 | 0.00 | 0.37 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.37 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.37 |  |
| 5000 | -1 | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  | 1 | 0.00 | 0.00 | 0.00 | 0.01 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.01 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.01 |  |
Medium | 200 | -1 | 0.00 | 0.00 | 0.00 | 0.04 | 0.47 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.04 | 0.47 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.04 | 0.47 | 0.00
|  | 0 | 0.00 | 0.00 | 0.00 | 0.12 | 0.06 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.12 | 0.06 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.12 | 0.06 | 0.00
|  | 1 | 0.00 | 0.00 | 0.00 | 0.63 | 0.51 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.63 | 0.51 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.63 | 0.51 | 0.00
| 1000 | -1 | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  | 0 | 0.00 | 0.00 | 0.00 | 0.01 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.01 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.01 |  |
|  | 1 | 0.00 | 0.00 | 0.00 | 0.10 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.10 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.10 |  |
| 5000 | -1 | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  | 0 | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  | 1 | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
High | 200 | -1 | 0.00 | 0.00 | 0.00 | 0.04 | 0.46 | 0.00