arXiv:1511.08404v2 [] 19 Aug 2016

Improved Precision in the Analysis of Randomized
Trials with Survival Outcomes, without Assuming
Proportional Hazards
Iván Dı́az1 , Elizabeth Colantuoni2 , Daniel F. Hanley3 , and Michael Rosenblum2
1

Division of Biostatistics and Epidemiology, Weill Cornell Medicine, New York, NY, USA
of Biostatistics, Johns Hopkins Bloomberg School of Public Health, Baltimore, MD, USA.
3 Division of Brain Injury Outcomes, Johns Hopkins Medical Institutions, Baltimore, MD, USA

2 Department

August 22, 2016
Abstract
We present a new estimator of the restricted mean survival time in randomized trials where
there is right censoring that may depend on treatment and baseline variables. The proposed
estimator leverages prognostic baseline variables to obtain equal or better asymptotic precision compared to traditional estimators. Under regularity conditions and random censoring
within strata of treatment and baseline variables, the proposed estimator has the following features: (i) it is interpretable under violations of the proportional hazards assumption; (ii) it is
consistent and at least as precise as the Kaplan-Meier estimator under independent censoring;
(iii) it remains consistent under violations of independent censoring (unlike the Kaplan-Meier
estimator) when either the censoring or survival distributions are estimated consistently; and
(iv) it achieves the nonparametric efficiency bound when both of these distributions are consistently estimated. We illustrate the performance of our method using simulations based on
resampling data from a completed, phase 3 randomized clinical trial of a new surgical treatment for stroke; the proposed estimator achieves a 12% gain in relative efficiency compared to
the Kaplan-Meier estimator. The proposed estimator has potential advantages over existing approaches for randomized trials with time-to-event outcomes, since existing methods either rely
on model assumptions that are untenable in many applications, or lack some of the efficiency
and consistency properties (i)-(iv). We focus on estimation of the restricted mean survival time,
but our methods may be adapted to estimate any treatment effect measure defined as a smooth
contrast between the survival curves for each study arm. We provide R code to implement the
estimator.

Keywords: Covariate adjustment; Efficiency; Targeted minimum loss based estimation; Random censoring.

1

1

Introduction

A standard approach to analyze clinical trials with survival outcomes is to estimate the survival
curve in each study arm using the Kaplan-Meier estimator. This approach assumes that the censoring time is independent of the event time for each study arm. This assumption would typically be
false, e.g., if high baseline disease severity is prognostic for both earlier drop-out and earlier time
to death. To accommodate this type of situation, we focus throughout on the weaker assumption
of random censoring, defined as censoring being independent of the event time within strata of the
study arm and baseline variables. This assumption allows informative censoring, i.e., censoring
correlated with the event time in a manner fully explained by study arm and baseline variables.
The Kaplan-Meier estimator is an unadjusted estimator, that is, it ignores baseline variables.
Unadjusted estimators have the following potential drawbacks: (a) they can yield inconsistent estimators of the survival function under informative censoring as discussed above, and (b) even under
independent censoring (defined as censoring independent of event time and baseline variables, for
each arm), they are inefficient if baseline variables are prognostic of the outcome. Appropriate adjustment for baseline variables provides an opportunity to avoid these drawbacks by (a) providing
consistent estimators of the treatment effect under random censoring and consistent estimation of
either the censoring or survival distribution, and (b) increasing efficiency (i.e., asymptotic precision) of the estimators, thereby decreasing the required sample size and saving resources.
A commonly used method to analyze randomized trials with survival outcomes is the proportional hazards model (Cox, 1972). A drawback of this approach is that the treatment effect estimate
becomes uninterpretable and may be misleading under violations of the proportional hazards assumption (Schemper, 1992; Tian et al., 2014). These violations are not necessarily easy to detect,
and can lead to false conclusions in an otherwise well designed and executed randomized trial.
Several alternatives are available to define the effect of assignment to treatment versus control
on a survival outcome. We focus on estimation of the marginal (i.e., unconditional) treatment
effect defined as the difference between the restricted mean survival time (RMST) in the two study
arms. The RMST is the expected survival time restricted to (i.e., truncated at) a time τ . This
parameter has a model-free, clinically meaningful interpretation (Chen and Tsiatis, 2001; Royston
and Parmar, 2011; Zhao et al., 2012; Tian et al., 2014). For example, if τ = 180 days, then a 14
day improvement in RMST due to treatment means 2 more weeks alive on average during the first
6 months; this may be more directly interpretable than a hazard ratio.
Although we focus on the RMST, our methods may be adapted to estimation of any smooth
contrast between the marginal survival curves for each study arm. For example, the difference
between the restricted median survival times may be of interest, e.g., for heavy-tailed survival
times. An alternative goal, not considered here, is to estimate a conditional effect of treatment, i.e.,
a contrast between distributions conditioned on the values of certain baseline variables. Though
our ultimate goal is to estimate unconditional treatment effects, we harness information in baseline
variables both to handle informative censoring (which could lead to bias if ignored) and to adjust
for chance imbalances between study arms (to improve precision).
Our methods assume the outcome is observed on a discrete time scale. If the outcome is
measured on a continuous time scale, our proposal may still be used by finely discretizing time. In
our motivating example, we observe time at the day level, and consider a period of τ = 180 days.
Existing estimators for our problem can be broken into the following three families: methods
based on modeling the survival function conditional on treatment and baseline variables (referred to
2

as outcome regression models), methods based on modeling the censoring probability conditional
on treatment and baseline variables, and methods called doubly robust estimators that combine
these models. We describe these in Section 2 where we show existing methods lack one or more
of the features (i)-(iv) in the abstract.
We propose a new doubly robust estimator of the RMST with all the properties in the abstract,
which is derived using the general targeted minimum loss based estimation (TMLE) framework
of van der Laan and Rubin (2006). Our estimator combines key ideas from Moore and van der
Laan (2009a), Rotnitzky et al. (2012), and Gruber and van der Laan (2012), as we describe in
Section 6. To the best of our knowledge, our proposed estimator is the first to achieve properties
(i)-(iv) simultaneously for our problem.
In Section 2, we review commonly used methods for our problem and highlight their strengths
and limitations. In Section 3, our motivating application is presented: the analysis of a completed
Phase III randomized trial of a new treatment for stroke. In Sections 4 and 5, we define our estimation problem and present estimators from related work, respectively. In Section 6, we present our
new estimator. Simulation studies are presented in Section 7, based on our motivating application.
These simulation studies demonstrate that our estimator can lead to substantial improvements in
efficiency. We conclude with a brief discussion and directions of future research.

2

Related Work

Various methods have been proposed that satisfy some but not all of the properties (i)-(iv) from the
abstract. Zhang (2014) propose a method to estimate the survival function at a single time point, by
using a linear working model for the outcome. A working model is defined as a statistical model
used to construct an estimator, but that is not necessarily assumed to be correctly specified. Unlike
our proposal, the method of Zhang (2014) requires censoring to be independent of the event time
within each study arm. In addition, unlike our estimator, their method requires outcome regression
models to be linear (and not, e.g., logistic). Parast et al. (2014) propose an estimator that adjusts for
covariates through a kernel regression of the outcome on a one-dimensional dimension reduction
defined as the linear predictor of a proportional hazards model. Their estimator achieves properties
(i)-(ii), but not (iii)-(iv).
Lu and Tsiatis (2011, Section 3) give a general approach for constructing estimators using
longitudinal data. Their approach, if it were applied to estimate the RMST, would have properties
(i)-(ii), but not (iii)-(iv). Unlike here, they require the censoring distribution to be known. An
advantage of their estimator over our proposal is the incorporation of post-baseline covariates to
adjust for time dependent confounding (i.e., when censoring and the event time have time-varying
common causes). Though we do not address this problem, the techniques in this paper may be
generalized to accommodate such scenarios. Stitelman et al. (2011) also handle time dependent
confounding for survival outcomes; their estimator has properties (i), (iii), (iv) but not (ii).
Methods based on estimating the censoring distribution (Cole and Hernán, 2004; Xie and Liu,
2005; Rotnitzky and Robins, 2005) are consistent under correct specification of the censoring
model. However, these estimators are typically not as efficient as the covariate adjusted estimators
described later, and thus do not fully leverage the often expensive data collected in a clinical trial.
Under consistent estimation of the censoring and outcome distributions at rate faster than n1/4 ,
doubly robust estimators are asymptotically efficient in the nonparametric model that only assumes
3

treatment is assigned independent of baseline variables (Robins and Rotnitzky, 1992; Hubbard
et al., 2000; van der Laan and Robins, 2003; Moore and van der Laan, 2009a; Stitelman et al.,
2011). Under outcome regression model misspecification but correct censoring model specification, doubly robust estimators remain consistent but they can be inefficient with variance larger
than the variance of inverse probability weighted estimators. Under independent censoring, they
can also have variance larger than unadjusted alternatives such as the Kaplan-Meier estimator. (A
large sample illustration, for a related problem, is in the Web Appendix of Dı́az et al., 2015). This
is problematic since there is no guarantee that the effort placed in constructing adjusted estimators
will lead to improved precision in estimation of the treatment effect.
The efficiency theory that we develop has roots in the work of Pfanzagl and Wefelmeyer (1985);
Robins and Rotnitzky (1992); Robins et al. (1994); Bickel et al. (1997); Hahn (1998); Scharfstein
et al. (1999); Bang and Robins (2005), among others who laid the foundation for locally efficient
estimation of causal effects. Their methods have been extended to incorporate enhanced efficiency
properties, e.g., by Tan (2006); van der Laan and Rubin (2006); Zhang et al. (2008); Tsiatis et al.
(2008); Cao et al. (2009); Tan (2010); Rotnitzky et al. (2012); Gruber and van der Laan (2012).
We use the general framework of van der Laan and Rubin (2006) and Gruber and van der
Laan (2012) to construct a targeted minimum loss based estimator of the RMST that satisfies the
properties described in the abstract. Targeted minimum loss based estimation (TMLE) of the effect
of treatment on binary, continuous, and time to event outcomes in randomized trials was discussed
in Moore and van der Laan (2009b) and Stitelman et al. (2011). Based on the work of Gruber and
van der Laan (2012), Dı́az et al. (2015) proposed a TMLE for ordinal outcomes with enhanced
efficiency properties analogous to those of the estimator we propose in this manuscript.
The general estimation approach of Moore and van der Laan (2009a) has properties (i), (iii),
(iv), but not (ii). Our main innovation is to enhance this approach so that it also achieves (ii).
The enhancement is not trivial to achieve, and relies on the general strategy from Rotnitzky et al.
(2012) and Gruber and van der Laan (2012), whose work builds on enhanced efficiency methods
described above.

3

Motivating Application: CLEAR III Trial

The CLEAR III trial (Clot Lysis: Evaluating Accelerated Resolution of Intraventricular Hemorrhage Phase III) is a completed Phase III multicenter randomized controlled trial. 500 individuals
with intraventricular hemorrhage (IVH) were randomized with equal allocation to receive either
alteplase (treatment) or saline (control) for IVH removal. The primary outcome was defined as a
score of 3 or less on the modified Rankin Scale (mRS) of functional disability at 180 days (where
smaller values correspond to better function). The proportions with mRS at most 3 at 180 days
were 48% vs. 45% in the treatment vs. control arms, respectively. A key secondary outcome was
all-cause mortality at 180 days: 18% and 29% of patients experienced death by 180 days in the
treatment vs. control arms, respectively. Though covariate adaptive randomization was used to
assign study arms, we ignore this in illustrating our method.
We reanalyzed data from the CLEAR III trial by defining the outcome as time to death in days
from randomization, and the treatment effect as the difference in the RMST (at τ = 180 days)
comparing treatment versus control arms. Figure 1 displays the estimated Kaplan-Meier survival
curve for each study arm. Using the Kaplan-Meier estimator, the unadjusted estimate of the differ4

ence in the RMST is 14.9 days (with standard error 5.5). The CLEAR III investigators identified
several baseline variables believed to be prognostic for mortality; these include age, the Glasgow
Coma Score (GCS), the National Institute of Health Stroke Scale (NIHSS) score, intracerebral
hemorrhage (ICH) location (thalamus vs. other) and ICH volume. After adjusting for these baseline variables using our proposed method described in Section 6, the estimated difference in the
RMST is 14.6 days (with standard error 5.2). Our adjusted estimator yields an estimated variance
that is roughly 12% smaller than the (unadjusted) RMST difference based on the Kaplan-Meier
estimator. In the context of the CLEAR III trial, such a precision gain would allow a reduction
by approximately 60 (out of original 500) patients in the required sample size to achieve a desired
power, if a Wald-test were used based on the adjusted versus the unadjusted estimator.

1.0

Proportion of Patients Alive

Treatment
Control

0.9

0.8

0.7

0.6
0

30

60

90

120

150

180

Time (in days)

Figure 1: Kaplan-Meier survival curves for the treatment arm (solid line) and control arm (dashed
line) of the CLEAR III trial

5

4
4.1

Data Structure, RMST Parameter, and Identification
Observed Data Structure for Each Participant

Assume K equally spaced time points t = {1, . . . , K}, e.g., representing days, at which participants are monitored. Let T denote a discrete, time-to-event outcome taking values in {1, . . . , K} ∪
{∞}, where T = ∞ represents no event occurring during times 1, . . . , K. Let C ∈ {0, . . . , K}
denote the censoring time defined as the time at which the participant is last observed in the study;
if a participant remains on study through time point K, we let C = K, which represents administrative censoring. Let A ∈ {0, 1} denote study arm assignment, and let W denote a vector of
baseline variables. The observed data vector for each participant is O = (W, A, ∆, T̃ ), where
T̃ = min(C, T ), and ∆ = 1{T ≤ C} is the indicator that the participant’s event time is observed
(uncensored). Here 1(X) is the indicator variable taking value 1 if X is true and 0 otherwise.
We assume the observed data vector for each participant i, denoted Oi = (Wi , Ai , ∆i , T̃i ), is an
independent, identically distributed draw from the unknown joint distribution P0 on (W, A, ∆, T̃ ).
We assume P0 ∈ M, where M is the nonparametric model defined as all continuous densities
on O with respect to a dominating measure ν such that A is independent of W , which holds by
randomization. Our asymptotic results are in the limit as sample size n goes to infinity, with the
number of time points K being fixed.
We can equivalently encode a single participant’s data vector O using the following longitudinal data structure:
O = (W, A, R0 , L1 , R1 , L2 . . . , RK−1 , LK ),
(1)
where Rt = 1{T̃ = t, ∆ = 0} and Lt = 1{T̃ = t, ∆ = 1}, for t ∈ {0, . . . , K}. The sequence
R0 , L1 , R1 , L2 . . . , RK−1 , LK in the above display consists of all 0’s until the first time that either
the event is observed or censoring occurs, i.e., time t = T̃ . In the former case Lt = 1; otherwise
Rt = 1. For a random variable X, we denote its history through time t as X̄t = (X0 , . . . , Xt ). For
a given scalar x, the expression X̄t = x denotes element-wise equality. The corresponding vector
(1) for participant i is denoted by (Wi , Ai , R0,i , L1,i , R1,i , L2,i . . . , RK−1,i , LK,i ).
Define the following indicator variables for each t ≥ 1:
It = 1{R̄t−1 = 0, L̄t−1 = 0},

Jt = 1{R̄t−1 = 0, L̄t = 0}.

The variable It is the indicator based on the data through time t − 1 that a participant is at
risk of the event being observed at time t; in other words, It = 1 means that all the variables
R0 , L1 , R1 , L2 ..., Lt−1 , Rt−1 in the data vector (1) equal 0, which makes it possible that Lt = 1.
Analogously, Jt is the indicator based on the outcome data through time t and censoring data
before time t that a participant is at risk of censoring at time t. By convention we let J0 = 1.
Define the hazard function for survival at time m ∈ {1, . . . , K}:
h(m, a, w) = P0 (Lm = 1|Im = 1, A = a, W = w),
among the population at risk at time m within strata of study arm and baseline variables. Similarly,
for the censoring variable C, define the censoring hazard at time m ∈ {0, . . . , K}:
gR (m, a, w) = P0 (Rm = 1|Jm = 1, A = a, W = w).
6

We use the notation gA (a, w) = P0 (A = a|W = w) and g = (gA , gR ). Let pW denote the
marginal distribution of the baseline variables W . We add the subscript 0 to pW , g, h to denote
the corresponding quantities under P0 . The joint distribution P0 on the observed data vector O =
(W, A, ∆, T̃ ) is completely characterized by the components p0 , g0 , h0 , i.e., P0 = (p0 , g0 , h0 ).

4.2

RMST Parameter Definition in Terms of Potential Outcomes

Define the potential outcomes Ta : a ∈ {0, 1} as the event times that would have been observed
had study arm assignment A = a and censoring time C = K been externally set with probability
one. For a restriction time τ ∈ {1, . . . , K} of interest, the target estimand is the difference between
the restricted mean survival time setting study arm to a = 1 versus a = 0:
θc = E{min(T1 , τ ) − min(T0 , τ )}.
The superscript c denotes a causal parameter, that is, a parameter of the distribution of the potential outcomes T1 and T0 . We prove in the Supplementary Materials that E{min(Ta , τ )} =
P
τ −1 c
c
t=0 S (t, a), where S (t, a) = P (Ta > t) is the survival probability corresponding to the potential outcome under assignment to arm A = a. As a result, θc may be expressed as
c

θ =

τ −1
X

{S c (t, 1) − S c (t, 0)},

(2)

t=1

since S c (0, a) = 1 for a ∈ {0, 1}.

4.3

Identification of RMST Parameter θc in Terms of Observed Data Generating Distribution P0

We show how the RMST parameter θc , which is defined above in terms of potential outcomes, can
be equivalently expressed as a function θ of the observed data distribution P0 (W, A, ∆, T̃ ), under
the assumptions (a)-(d) below. This is useful since the potential outcomes are not always observed,
in contrast to the observed data vector (W, A, ∆, T̃ ) for each participant (whose distribution we can
make direct statistical inferences about); we refer to θ as a statistical parameter, which is shorthand
for saying it is a mapping from the observed data distribution P ∈ M to R.
Define the following assumptions:
(a) T = 1(A = 0)T0 + 1(A = 1)T1 (consistency);
(b) A is independent of (Ta , W ), for each a ∈ {0, 1} (randomization);
(c) C is independent of Ta conditional on (A, W ), for each a ∈ {0, 1} (random censoring);
(d) gA,0 (a, w) > 0 and gR,0 (t, a, w) < 1 whenever the P0 -density of W is positive at W = w, for
each a ∈ {0, 1} and t ∈ {0, . . . , τ − 1} (positivity assumption).
We make assumptions (a)-(d) throughout the manuscript. Assumption (a) connects the potential
outcomes to the observed outcome. Assumption (b) holds by design in a randomized trial. Assumption (c), which is similar to that in Rubin (1987), means that censoring is random within
7

strata of treatment and baseline variables (which we abbreviate as “random censoring”). Assumption (d) states that each treatment arm has a positive probability, and that every time point has a
hazard of censoring smaller than one, within each baseline variable stratum W = w with positive
density under P0 .
Denote the survival function for T at time t ∈ {1, . . . , τ − 1} conditioned on study arm a and
baseline variables w by
S(t, a, w) = P (T > t|A = a, W = w).
(3)
Similarly, define the following function of the censoring distribution:
G(t, a, w) = P (C ≥ t|A = a, W = w).

(4)

Under assumptions (a)-(d), we have T ⊥⊥C|A, W and therefore S(t, a, w) and G(t, a, w) have the
following product formula representations:
S(t, a, w) =

t
Y

{1 − h(m, a, w)};

G(t, a, w) =

t−1
Y

{1 − gR (m, a, w)}.

(5)

m=0

m=1

The potential outcome survival function S c (t, a) can be equivalently represented in terms of
the observed data distribution as
S(t, a) = EpW

t
Y

{1 − h(m, a, W )},

(6)

m=1

for t ∈ {1, . . . , K}, a ∈ {0, 1}; equality of S c (t, a) and the above display follows from (5) and
S c (t, a) = P (Ta > t) = EpW P (Ta > t|W ) = EpW P (T > t|A = a, W ) = EPW S(t, a, W ),
where the third equality above follows from (a) and (b).
It follows from (2) that the causal parameter θc is equal to the following statistical parameter:
( t
)#
( t
)
"
τ −1
Y
Y
X
E PW
{1 − h(m, 1, W )} − EPW
{1 − h(m, 0, W )}
.
(7)
θ=
t=1

m=1

m=1

Our goal is to estimate θ based on n independent, identically distributed observations Oi =
(Wi , Ai , ∆i , T̃i ) drawn from P0 . We construct an estimator with properties (i)-(iv) in the abstract.
Since the parameter of interest θ is defined as a function of (pW , h) through (7), a natural estimation strategy would be to plug estimates of pW and h in these formulas. Estimators constructed
in this way are called substitution or plug-in estimators, and have the advantage that they remain
within bounds of the parameter space; this is desirable in estimation of probabilities and other
bounded parameters such as the RMST. Our proposed estimator is a substitution estimator.
Define independent censoring to be C⊥⊥(Ta , W ) | A, for each a ∈ {0, 1}. This is a stronger
(more restrictive) assumption than random censoring. We refer to a censoring mechanism G(t, a, w)
as non-informative if it does not depend on w, and as informative if it depends on w. Under independent censoring, G(t, a, w) is non-informative.

8

5

Several Estimators of RMST (θ) from Related Work

5.1

Unadjusted Estimators of θ: Kaplan-Meier and Inverse Probability Weighted

ThroughoutQthis subsection only, we additionally assume independent censoring. This implies
S(t, a) = tm=1 {1 − h(m, a)} for h(m, a) = P (Lm = 1|Im = 1, A = a), as proved in the
Supplementary Material. Therefore, we have the following simpler representation of θ:
" t
#
τ −1 Y
t
X
Y
θ=
{1 − h(m, 1)} −
{1 − h(m, 0)} .
(8)
t=1

m=1

m=1

The Kaplan-Meier estimator for S(t, a) is defined as
Ŝkm (t, a) =

t 
Y

Pn
1−

m=1

1{Lm,i = 1, Im,i = 1, Ai =
i=1P
n
i=1 1{Im,i = 1, Ai = a}

a}


,

(9)

where we set the above fraction to be 0 if the denominator is 0. The right side of (9) was obtained
by substituting the empirical counterpart of each h(m, a) in the formula for S(t, a) above. The
corresponding Kaplan-Meier estimator of θ is defined analogously, as
τ −1
X
{Ŝkm (t, 1) − Ŝkm (t, 0)}.
θ̂km =
t=1

Since θ̂km is a smooth function of at most 4(τ − 1) empirical means, the delta method (Theorem
3.1 van der Vaart, 1998) implies
√

n

1 X
n(θ̂km − θ) = √
D (Oi ) + oP (1),
n i=1 km

where, for an observation O,

τ −1 X
t 
X
(2A − 1)Im S(t, A)
Dkm (O) = −
{Lm − h(m, A)}
g
(A)G(m,
A)
S(m,
A)
A
t=1 m=1

(10)

is the influence function of θ̂km . The above influence function Dkm may be derived from Lemma 1
below, noting that the Kaplan-Meier estimator is the maximum likelihood estimator in the nonparametric model where only (A, ∆, T̃ ) (and not W ) are measured, and therefore it is asymptotically
linear with influence√function equal to the efficient influence function in the model for (A, ∆, T̃ ).
As a consequence, n(θ̂km − θ) converges to a mean zero normal distribution with asymptotic
variance Var(Dkm (O)).
This estimator is consistent and efficient for the case where only (A, ∆, T̃ ) is observed. However, the unadjusted estimator is generally inefficient for the case where (W, A, ∆, T̃ ) is observed.
Intuitively, this is because the unadjusted estimator fails to leverage the prognostic information in
baseline variables W . Furthermore, under the less restrictive random censoring assumption, the

9

unadjusted estimator will generally not be consistent, while adjusted estimators remain consistent
if the censoring distribution is consistently estimated.
We next define the unadjusted, inverse probability weighted (IPW) estimator.
Let Ĝkm (t, a)
Q
denote the Kaplan-Meier estimator of the censoring distribution G(t, a) = t−1
{1
− gR (m, a)},
m=0
defined as
Pn

t−1 
Y
1{Rm,i = 1, Jm,i = 1, Ai = a}
i=1P
.
Ĝkm (t, a) =
1−
n
i=1 1{Jm,i = 1, Ai = a}
m=0
The unadjusted IPW estimator of S(t, a) is defined as
n

Ŝipw (t, a) =

1 X 1{Ai = a, R̄t−1,i = 0, L̄t,i = 0}
,
n i=1
ĝA (a)Ĝkm (t, a)

where ĝA (a) denotes the sample mean of 1{A = a}. Zhao et al. (2016) show that
√
n{Ŝipw (t, a) − Ŝkm (t, a)} = oP (1),
√
i.e., the asymptotic distributions
Pτ −1of these two estimators are equal up to oP (1/ n), and the estimator of θ given by θ̂ipw = t=1 {Ŝipw (t, 1) − Ŝipw (t, 0)} also satisfies
n
√
1 X
D (Oi ) + oP (1).
n(θ̂ipw − θ) = √
n i=1 km

√
√
An important consequence is that both n(θ̂km − θ) and n(θ̂ipw − θ) converge in distribution to
2
2
), where σkm
= Var(Dkm (O)).
N (0, σkm

5.2

Covariate Adjustment: Inverse Probability Weighted Estimators

Consider an estimator ĥ(t, a, w) for the outcome hazard function h(t, a, w), for example, based on
fitting a Cox proportional hazards model that conditions on A = a and W = w, and that uses
the Nelson-Aalen estimator for the baseline hazard. If the proportional hazards model is correct,
the corresponding substitution estimator based on (7) is consistent for θ. However, under model
misspecification, this substitution estimator will generally be inconsistent. This is particularly
problematic for randomized trials with independent censoring, in which an unadjusted, consistent
estimator can be obtained through the Kaplan-Meier survival function.
As an alternative, an adjusted, inverse probability weighted (IPW) estimator of S(t, a) (where
here and below we use the definition of S(t, a) in (6)) is given by
n

Ŝadj,ipw (t, a) =

1 X 1{Ai = a, R̄t−1,i = 0, L̄t,i = 0}
,
n i=1
ĝA (a, Wi )Ĝ(t, a, Wi )

where ĝA (a, w) and Ĝ(t, a, w) are estimators of gA (a, w) and G(t, a, w), respectively.
The adjusted
Pτ −1
IPW estimator of the τ -restricted mean survival time, given by θ̂adj,ipw = t=1 {Ŝadj,ipw (t, 1) −
Ŝadj,ipw (t, 0)}, is consistent for θ if the estimators (ĝA , ĝR , ĥ) are consistent at a fast enough rate.

10

We focus on the case where ĝR (t, a, w) and ĝ(a, w) are estimators having saturated terms for
time, treatment, and their interaction. We define such estimators as estimators that satisfy
Pn

i=1

Jm,i {Rm,i − ĝR (m, Ai , Wi )} =

Pn

i=1

Ai Jm,i {Rm,i − ĝR (m, Ai , Wi )} =

Pn

i=1 {Ai

− ĝA (1, Wi )} = 0,

for m = 0, . . . , K − 1. An example of estimators with saturated terms for time, treatment and their
interaction is given by maximum likelihood estimators for logistic regression models of gR (t, a, w)
and gA (a, w), respectively, that include at least an intercept, main terms for A and each time t ∈
{1, . . . , K}, and interaction terms for A by each time t ∈ {1, . . . , K}; arbitrary and data-adaptive
additional terms involving t, a, w can be included in these models.
When using Ĝ(t, a, w) and ĝ(a, w) estimated with such saturated terms, under independent
censoring, Ŝadj,ipw (t, a) is consistent for S(t, a). If these models contain no terms involving w,
then Ŝadj,ipw (t, a) and Ŝipw (t, a) (both using the corresponding estimators Ĝ and ĝ) are identical. Williamson et al. (2014) use this observation to show that, under independent censoring, the
asymptotic variance of Ŝadj,ipw (t, a) is smaller or equal than the asymptotic variance of Ŝipw (t, a);
this, together with the delta method, implies
√
2
n(θ̂adj,ipw − θ) → N (0, σadj,ipw
),
2
2
where σadj,ipw
≤ σkm
. (Throughout, → indicates convergence in distribution as n goes to infinity.)
A weakness of θ̂adj,ipw is that under informative censoring, it will generally be inconsistent if the
model for G is misspecified. This motivates considering double robust estimators described in
Section 5.3.
Adjusted estimators (such as the adjusted IPW above) often involve fitting a parametric model
for gA , even though gA is known by design in a randomized trial. Intuitively, the purpose of this
model fit is to capture chance imbalances of the baseline variables W between study arms for a
given data set; these imbalances can then be adjusted to improve efficiency. The general theory
underlying efficiency improvements through estimation of known nuisance parameters such as gA
is presented, e.g., by Robins et al. (1994) and van der Laan and Robins (2003).

5.3

Augmented Inverse Probability Weighted (Double Robust) Estimator

We start by presenting the efficient influence function for estimation of θ in model M. The following lemma may be proved by applying the delta method to the definition of θ in (7) and using
the efficient influence function of S(t, a) as presented in Moore and van der Laan (2009a):
Lemma 1. The efficient influence function for estimating θ in the model M is
D(O) =

τ −1
X


Im Z(m, A, W ) {Lm − h(m, A, W )} + S(m, 1, W ) − S(m, 0, W ) − θ,

(11)

m=1

where Z(m, A, W ) = Z1 (m, A, W ) − Z0 (m, A, W ), and
Za (m, A, W ) = −

τ −1
X

S(t, a, W )
1{A = a}
.
g (a, W )G(m, a, W ) S(m, a, W )
t=m A
11

(12)

For conciseness, we suppress the dependence of D on g and h. The function D has two
important properties for estimation of θ. First, it is a doubly robust estimating function, i.e., for
given estimators ĥ and ĝ of h and g = (gA , gR ), respectively, the estimator formed by solving for θ
in the following estimating equation is consistent if at least one of h or g is estimated consistently
(while the other converges to a limit that may be incorrect):
0=

n X
τ −1  X
τ −1
X
i=1

m=1

t=m

−(2Ai − 1)Im,i

Ŝ(t, Ai , Wi )

ĝA (Ai , Wi )Ĝ(m, Ai , Wi ) Ŝ(m, Ai , Wi )

{Lm,i − ĥ(m, Ai , Wi )}



+ Ŝ(m, 1, Wi ) − Ŝ(m, 0, Wi ) − θ , (13)
Q
Qm−1
0
0
where Ŝ(m, a, w) = m
m0 =1 {1 − ĥ(m , a, w)} and Ĝ(m, a, w) =
m0 =0 {1 − ĝR (m , a, w)}. This
estimator is often referred to as the augmented IPW estimator, and we denote it by θ̂aipw . The
double robustness property is desirable since it guarantees that improper adjustment for covariates
through a misspecified working model for h still leads to a consistent estimator of θ in randomized
trials with random censoring if gR is consistently estimated.
Second, the efficient influence function (11) characterizes the information bound for estimation
of θ in the model M (Bickel et al., 1997). Specifically, under consistent estimation of h and g at
a sufficiently fast rate, θ̂aipw has variance smaller or equal to that of any regular, asymptotically
linear estimator of θ in M (a property called local efficiency); in this case, if independent censoring
holds, then θ̂aipw has equal or smaller asymptotic variance compared to θ̂km . Unfortunately, under
misspecification of the model for h, the estimator of h will generally be inconsistent, which could
lead to θ̂aipw having worse asymptotic efficiency than θ̂km under independent censoring. In other
words, the added robustness of θ̂aipw may come at the price of lower efficiency compared to θ̂km .
This motivates the question of whether this added robustness can be achieved at no cost. I.e., can
we construct a doubly robust estimator with the added guarantee of equal or better asymptotic
precision as θ̂km if independent censoring holds? We construct such an estimator below.

6

Proposed Estimator

We develop a doubly robust estimator, denoted θ̂adj,eff that has properties (i)-(iv) of the abstract.
We will sometimes use the following modified representation of the data set:
{(m, Wi , Ai , Jm,i , Rm,i , Im+1,i , Lm+1,i ) : m = 0, . . . , K − 1; i = 1, . . . , n}.

(14)

This data set is referred to as the long form, and the original data set
{(Wi , Ai , ∆i , T̃i ) : i = 1, . . . , n}

(15)

is referred to as the short form. An observation in the long form data set is a vector of the form
(m, Wi , Ai , Jm,i , Rm,i , Im+1,i , Lm+1,i ). We define the following auxiliary variables:

12

M (W ) =

τ −1 
X
S(t, 1, W )
t=1


S(t, 0, W )
+
,
gA (1, W )
gA (0, W )

τ −1
X
(2A − 1) S(t, A, W )
1
H(m, A, W ) = −
.
g
(A,
W
)
S(m,
A,
W
)
G(m
+
1,
A,
W
)
A
t=m+1

(16)
(17)

These auxiliary covariates are constructed from a decomposition of the efficient influence function
targeting improved efficiency of the estimators as described after Theorem 2 below.
Our proposed estimator θ̂adj,eff is defined as the output of the following algorithm (which follows the TMLE template):
Step 1. Initial estimators. Obtain initial estimators ĝA , ĝR , and ĥ of gA , gR , and h, respectively.
Let p̂W denote the empirical distribution of W .
Step 2. Iteratively update estimates of h and g. Initialize l = 0, and let ĥl = ĥ, ĝAl = ĝA , ĝRl = ĝR .
(a) Update ĥl . Let Ŝ l , Ĝl denote S, G after substituting ĥl , ĝRl , p̂W for h, gR , pW in (5)
and (7). Augment each observation in the long form data set (14) by two covariates
Ẑal (m, Ai , Wi ) : a ∈ {0, 1}, where each Ẑal (m, Ai , Wi ) is constructed by substituting
the estimates Ŝ l , Ĝl , ĝAl evaluated at (m, Ai , Wi ) in (12). Estimate the parameter vector
 = (1 , 0 ) in the logistic hazard submodel hl for h:
logit hl (m, a, w) = logit ĥl (m, a, w) + 1 Ẑ1l (m, a, w) + 0 Ẑ0l (m, a, w),

(18)

by computing the following maximum likelihood estimator:
ˆ = arg max


n X
τ −1
X

o
n
Im,i log hl (m, Ai , Wi )Lm,i (1 − hl (m, Ai , Wi ))1−Lm,i .

(19)

i=1 m=1

The maximizer ˆ can be computed using standard statistical software by a logistic
regression of Lm,i on the variables Z1l (m, Ai , Wi ), Ẑ0l (m, Ai , Wi ) among observations
with Im,i = 1 and m < τ in the long form data set (14), and using logit ĥl (m, Ai , Wi )
as an offset. Define ĥl+1 = hlˆ.
(b) Update ĝRl . Let Ĥ l (m, A, W ) denote H(m, A, W ) with Ŝ l , Ĝl , ĝAl substituted for
S, G, gA , respectively, in (17). Augment each observation in the long form data set
(14) by Ĥ l (m, Ai , Wi ). In the long form data set, estimate the parameter γ in the
following logistic regression submodel for gR (m, a, w):
l
logit gR,γ
(m, a, w) = logit ĝRl (m, a, w) + γ Ĥ l (m, a, w),

(20)

by logistic regression of Rm,i on the single covariate Ĥ l (m, Ai , Wi ) and with offset
logit ĝRl (m, A, W ) among observations with m < τ − 1 and Jm,i = 1. Denote the
corresponding maximum likelihood estimate of γ by γ̂.
13

(c) Update ĝAl . Let M̂ l (W ) denote M (W ) with ĝAl , Ŝ l substituted for the corresponding
components in (16). Augment each observation in the short form data set by M̂ l (Wi ).
In the short form data set, estimate the parameter ν in the following logistic regression
submodel for gA (1|w):
l
logit gA,ν
(1|w) = logit ĝAl (1|w) + ν M̂ l (w),

(21)

by logistic regression of Ai on the covariate M̂ l (Wi ) and with offset logit ĝAl (1|Wi )
among all participants i = 1, . . . , n. Denote the corresponding maximum likelihood
estimate of ν by ν̂.
l
l
Define ĥl+1 = hlˆ, ĝRl+1 = gR,γ̂
, and ĝAl+1 = gA,ν̂
.

Step 3. Update l = l + 1 and iterate the previous step until convergence. We stop at the first
iteration for which the sample mean of the squared differences of predictions based on ĥl ,
ĝRl , ĝAl between step l and step l + 1 is smaller or equal to 10−4 /n.
Denote ĥ? , ĝA? , and ĝR? the estimators obtained in the last iteration of the above algorithm, and
define the enhanced efficiency TMLE estimator of θ as
#
" n t
n Y
t n
τ −1
n
o
o
X
Y
X
X
1
1
1 − ĥ? (m, 1, Wi ) −
1 − ĥ? (m, 0, Wi ) .
(22)
θ̂adj,eff =
n i=1 m=1
n i=1 m=1
t=1
The above display is the substitution estimator of θ based on (7) where the sample means over baseline variables Wi , i = 1, . . . , n correspond to expectation with respect to the empirical distribution
of W .
We next consider ways to construct the initial estimators of h and g for step 1 above. The
outcome hazard function h may be estimated by running a prediction algorithm for the probability
of Lm = 1 as a function of A, W , and m among observations with Im = 1 in the long form data
set (14). The censoring hazard gR may be estimated by running an analogous prediction algorithm
of the probability that Rm = 1 as a function of A, W , and m among observations with Jm = 1.
The treatment mechanism gA may be estimated by fitting a parametric model for the probability
of A = 1 as a function of W in the short form data set. In a randomized trial, gA is set by design.
However, efficiency of the TMLE can be improved by estimating gA using, e.g., the proportion of
individuals in the treatment group, or a logistic regression model that contains baseline variables
and an intercept term.
Having defined the estimation algorithm, we now present our main results giving conditions
under which θ̂adj,eff is guaranteed to be at least as efficient as θ̂adj,ipw and θ̂km .
Theorem 1 (Equal or greater asymptotic efficiency compared to adjusted IPW estimator). Assume
(a)-(d), ĝA? , and ĝR? are n1/2 -consistent in L2 (P0 ) norm, and ĥ? converges to some limit h1 in L2 (P0 )
norm as n → ∞. Then we have the following convergence in distribution results:
√
√
2
2
n(θ̂adj,eff − θ) → N (0, σadj,eff
),
n(θ̂adj,ipw − θ) → N (0, σadj,ipw
),
2
2
where σadj,eff
≤ σadj,ipw
. In addition, if h1 equals the true h0 , then θ̂adj,eff achieves the semiparametric efficiency bound in M.

14

The consistency rates required in our
restrictive than necessary to
√ previous theorem are more
2
obtain the convergence in distribution n(θ̂adj,eff − θ) → N (0, σadj,eff ). In the Supplementary
Materials we present a more general theorem, along with its proof, which shows that this convergence holds under standard, less restrictive, doubly robust convergence assumptions on (ĝ ? , ĥ? ).
The following is our main result:
Theorem 2 (Equal or greater asymptotic efficiency compared to Kaplan-Meier estimator, under
independent censoring). Assume (a)-(d), independent censoring, and that ĥ? converges as in the
first sentence of Theorem 1. Assume also that ĝA? , ĝR? are estimated using models with saturated
terms for time, treatment, and their interaction. Then
√
√
2
2
),
),
n(θ̂adj,eff − θ) → N (0, σadj,eff
n(θ̂ipw − θ) → N (0, σkm
√
√
2
2
),
),
n(θ̂adj,ipw − θ) → N (0, σadj,ipw
n(θ̂km − θ) → N (0, σkm
2
2
2
where σadj,eff
≤ σadj,ipw
≤ σkm
.

These results guarantee that the proposed TMLE has asymptotic variance that never exceeds
that of the Kaplan-Meier estimator, under independent censoring (when the latter is consistent).
To the best of our knowledge, this is the first estimator to achieve the properties in the previous
theorems for our problem.
The algorithm above (i.e., steps 1-3 and the formula (22)) that generates the TMLE θ̂adj,eff
combines key ideas from Moore and van der Laan (2009a), Rotnitzky et al. (2012), and Gruber and
van der Laan (2012). Moore and van der Laan (2009a) present a TMLE algorithm for estimating
the survival difference S(t, 1) − S(t, 0) in the model M, which involves a step similar to 2a above;
this estimator has properties (i), (iii), and (iv) but not (ii). The crux of our approach to achieve
property (ii) (without sacrificing the other properties) is to augment the censoring and treatment
models through steps 2b and 2c. These augmented models use the covariates (16)-(17) that were
specifically constructed to achieve property (ii).
The idea of augmenting censoring and treatment models to achieve enhanced efficiency properties can be traced back at least to Robins et al. (1994). More recently, Rotnitzky et al. (2012) built
on this idea to construct an estimator that has equal or better asymptotic precision than a certain
parametric family of estimators including the adjusted IPW estimator, and then Gruber and van der
Laan (2012) showed how to do the same in the TMLE framework.
It is not trivial to determine precisely how to augment the censoring and treatment models
in order to guarantee (ii) holds. We explain the intuition for how we achieved this, which uses
general ideas from the above related work. Assume the conditions in Theorem 1. First, consider
the simpler case where the censoring and treatment distributions gR , gA are known. Define the
simplified TMLE to be as in steps 1-3 above, except omitting steps 2b and 2c and using the known
gR , gA in step 2a. The simplified TMLE’s influence function equals the influence function of the
adjusted IPW estimator minus the following expression (derived in the Supplementary Materials):
M (W ) [A − gA (1, W )] +

τ −2
X

Jm H(m, A, W ) [Rm − gR (m, A, W )] ,

m=0

for M and H the auxiliary variables in (16) and (17), respectively.
15

(23)

Second, consider the case where the censoring and treatment distributions gR , gA are unknown
and estimates ĝR , ĝA are used by the TMLE that involves all of steps 1-3, called the enhanced
efficiency TMLE (i.e., θ̂adj,eff ). The influence function for θ̂adj,eff equals that of the simplified
TMLE minus the latter’s projection on the tangent space Tg∗ spanned by the scores of the models
used to estimate gR , gA . (See van der Vaart (1998, Section 25.3) for background on tangent spaces
and projections as used here.) Subtracting off such a projection is helpful since it can only decrease
or leave unchanged the influence function’s variance, which equals the asymptotic variance of
the estimator. By augmenting the model for gR by H as in (20) in step 2b, the corresponding
score is the second term in (23); by augmenting the model for gA by M as in (21) in step 2c, the
corresponding score is the first term in (23). This implies (23) is in the tangent space Tg∗ , and
therefore the influence function for θ̂adj,eff is orthogonal to (23). Combining the above argument
with the last line in the previous paragraph, it follows that the influence function for θ̂adj,eff equals
the adjusted IPW influence function minus the latter’s projection on Tg∗ . Therefore, θ̂adj,eff has
asymptotic variance at most that of the adjusted IPW estimator, which gives the main conclusion
of Theorem 1. Theorem 2 then follows from the asymptotic equivalence of the adjusted IPW and
Kaplan-Meier estimators under the added assumptions of independent censoring and the model for
gR being saturated as described above. A more detailed argument that fleshes out and justifies the
above outline is in the Supplementary Materials.
In addition to the efficiency properties stated in the above results, our estimator inherits the
doubly robust property of the Moore and van der Laan (2009a) TMLE (i.e., the TMLE above
but without steps 2b and 2c). Under random censoring, this means that our proposal has two
opportunities to achieve consistency in estimating the causal effect, in contrast to the proportional
hazard model which relies exclusively on the assumption that the outcome regression is correctly
specified.
Under the assumptions of Theorem 1, for a consistent estimate σ̂adj,eff , a Wald-type confidence
√
interval θ̂adj,eff ± zα/2 σ̂adj,eff / n is guaranteed to have 1 − α asymptotic coverage probability.
If the initial estimators ĥ and ĝR in step 1 are M -estimators (e.g., if they are estimated through
maximum likelihood in parametric working models), the nonparametric bootstrap (Efron et al.,
1979) may be used to obtain a consistent estimate σ̂adj,eff (see Corollary 3.1 and 3.2 in Wellner
and Zhan, 1996). The performance of the nonparametric bootstrap is unknown when ĥ or ĝR are
data-adaptive estimators (e.g., if ĥ involves variable selection). The development of a consistent
variance estimator in this case remains an open question.

7
7.1

Simulation Study
Data Generating Distributions

Our simulated distributions are based on resampling data from the CLEAR III trial (described in
Section 3) in order to mimic key features of the trial. The baseline variables for each participant are
W = (W1 , W2 , W3 , W4 , W5 ) = (age, GCS, NIHSS score, ICH location, ICH volume), each scaled
to have mean 0 and standard deviation 1. The outcome is time until death in days, denoted by T .
Define the treatment effect θ to be the difference in the RMST at restriction time τ = 180 days.
We compare the performance of the following estimators defined above: Kaplan-Meier, unadjusted
16

IPW, adjusted IPW, augmented IPW (AIPW), and proposed TMLE. The adjusted IPW is referred
to simply as IPW. We ran simulations at two sample sizes: n = 500 and n = 2000.
We consider twelve data generating distributions, where we vary the correlation between T
and W (using 3 scenarios labeled A, B, C defined below), the treatment effect θ, and the censoring
mechanism. Some of these are set to mimic features of the CLEAR III data, as described below.
For each data generating distribution and sample size, we generated 10,000 simulated data sets,
and report on the empirical distribution of each estimator.
The key feature that determines the magnitude of precision gains for the adjusted estimators is
how correlated the baseline variables are with the outcome. To give a rough sense of the observed
correlations in the CLEAR III trial data set, we fit a logistic regression model for the hazard of death
that includes time, treatment, time by treatment interactions, and main effects for each baseline
variable in W = (W1 , W2 , W3 , W4 , W5 ). We report exp(β̂j ) : j = 1, 2, 3, 4, 5, where β̂j is the
estimated coefficient for Wj in the model fit. These were 1.5, 0.9, 1.7, 1.0, and 1.1, respectively.
We do not assume the above model is correct; we only used it to roughly measure correlations in
the CLEAR III data set. In our simulations, the data generating distributions for scenario A mimic
the above correlations; in scenarios B and C, the correlations are reduced as described below.
In each simulated trial, the data vector (W, A, T, C) for each participant is generated as an independent, identically distributed draw from a joint distribution P that satisfies assumptions (a)-(d).
This distribution will depend on the scenario A-C, censoring mechanism, and treatment effect θ.
For scenario A, each simulated participant’s data vector is generated by first resampling a participant with replacement from the CLEAR III data among the 491 patients who did not drop out,
and recording only his/her pair (W, T ). The resulting distribution preserves the empirical correlation between W and T from the CLEAR III trial. Such resampling may lead to more realistically
complex distributions than drawing from a regression model fit for T given W . In each scenario B
and C, we generate an initial (W, T ) as just described, and with probability 0.5 and 1, respectively,
we replace W by an independent draw (with replacement) from the marginal distribution of W in
the CLEAR III trial. The impact is that the correlation between T and W decreases as one proceeds from scenario A to scenario C, with scenario C having T and W independent (i.e., baseline
variables not prognostic for the outcome).
Next, we assign A independent of (W, T ) by a Bernoulli draw with probability 0.5 of being 1 or
0. This ensures that the randomization assumption (b) holds, and induces a distribution with θ = 0
(no treatment effect). We also consider distributions with a positive treatment effect (θ = 14.9);
this value was selected since it is the unadjusted Kaplan-Meier estimate of θ from the CLEAR III
data. To obtain such distributions, we generate (W, A, T ) as above except now if A = 1 we add
to T an independently generated draw from a χ2 distribution with mean µ = 56, where µ was
calibrated to achieve θ = 14.9.
For each scenario and treatment effect θ ∈ {0, 14.9}, given (W, A, T ) we generate C based
on either a non-informative or informative censoring model. Specifically, the censoring time C is
generated based on drawing from the corresponding distribution gR (t, a, w) given below, for each
t ∈ {0, . . . , T − 1} in turn:
Non-Informative gR : logit P [Rt = 1|Jt = 1, A, W ] = −5.5 + 0.007t;
Informative gR : logit P [Rt = 1|Jt = 1, A, W ] = −6.5 + 0.007t + 0.6W3 A + 0.3(W1 + W5 ).

17

The percentage of censored patients due to drop-out depends on the censoring model, treatment
effect and scenario. The non-informative and informative censoring distributions yield on average
62%–68% and 32%–38% censored patients due to drop-out, respectively (with values varying
within these intervals depending on θ and the scenario A-C). Though these are higher drop-out
rates than typically expected in a randomized trial, we evaluate performance under this substantial
censoring to illustrate that the estimators can have good performance even in this challenging
case. Assumptions (a)-(d) from Section 4.3 hold for all of our data generating distributions, while
independent censoring holds only under the non-informative censoring mechanism.

7.2

Estimators ĥ, ĝA , ĝR Used by the Adjusted Estimators

Each of ĥ, ĝA , ĝR is based on a logistic regression working model fit. The logistic regression
model for gA includes an intercept and a main term for each component of W . Since treatment
is assigned with P (A = 1|W ) = 0.5 for all simulated studies, the model for gA is correctly
specified. To account for censoring due to patient drop-out, the model for gR includes saturated
terms for time, treatment, and their interaction, in addition to main terms for W1 and W5 and a
treatment by W3 interaction. Therefore, the model for gR satisfies the condition in Theorem 2
and is correctly specified when the data is generated under either non-informative or informative
censoring. We use the long form data set to fit the model for gR , and set ĝR (m, a, w) = 0 for m
when the corresponding risk set is empty. The model for h consists of an intercept, main terms for
time (as a real value rather than categorical) and treatment, and treatment-time interaction, as well
as main terms for each component of W . For scenarios A and B, the data generating distribution
is based on resampling pairs W, T from the CLEAR III trial; therefore, the parametric model for h
is likely to be misspecified, which could easily occur in practice. However, in scenario C where T
is independent of W , the model for h is correctly specified.

7.3

Simulation Results

Tables 1 and 2 summarize the main results of our simulations, at sample sizes n = 500 and
n = 2000, respectively. Under non-informative censoring, where all the estimators are consistent,
all have relatively small bias. Under informative censoring and scenarios A and B (where the baseline covariates are prognostic for the outcome), the Kaplan-Meier and unadjusted IPW estimators
are more biased compared to their adjusted counterparts. The informative censoring distribution
depends on baseline variables (W1 , W3 , W5 ), which under scenarios A and B are also correlated
with the outcome, therefore causing confounding and non-negligible bias even at the large sample
size (n = 2000); the bias is relatively small for scenario C, since censoring and the outcome do not
share any common causes.
We measure relative efficiency as the ratio of mean squared error (RMSE) comparing the
Kaplan-Meier estimator to each of the other estimators. First consider scenario A, where the correlation between W and T mimics that from the CLEAR III trial data. The proposed estimator θ̂adj,eff
has relative efficiency gains in the range 12%–14% compared to the Kaplan-Meier estimator, under non-informative censoring. This means that a prespecified analysis plan using our proposed
TMLE could have required roughly 11% ≈ 1 − (1/1.12) fewer patients (55 out of 500) to achieve
the same power; alternatively, the trial sample size could be conservatively planned assuming no

18

precision gain, and then using θ̂adj,eff would lead to increased power if baseline variables are prognostic for the outcome. The efficiency gains of θ̂adj,eff are similar under both types of censoring
at sample size n = 500. However, the gains are larger under informative censoring at n = 2000
(reaching 53% in one case); this is due to the bias of the unadjusted estimators (which, unlike the
adjusted estimators, are inconsistent) making a substantial contribution to the mean squared error
at the larger sample size.
Next consider scenario C, where baseline variables are independent of T . For this scenario,
all of the estimators are consistent under both censoring distributions. There are small losses
in relative MSE for the adjusted estimators, due to their introducing unnecessary variability by
adjusting for baseline variables unrelated to the outcome. These precision losses are generally
smaller at the larger sample size n = 2000 (Table 2) compared to n = 500.
The unadjusted IPW estimator has similar bias and variance as the Kaplan-Meier estimator, as
predicted by theory. The adjusted IPW estimator is not as efficient as the double robust estimators
AIPW and TMLE in scenarios A and B.
The TMLE θ̂adj,eff has similar bias compared to the AIPW estimator. The TMLE variance is
slightly smaller, which translates into a relative MSE reduction of about 1% comparing the TMLE
to the AIPW estimator in scenarios A and B. The main advantages of the TMLE compared to
the AIPW estimator are: the TMLE has property (ii), i.e., guaranteed equal or better asymptotic
efficiency compared to the Kaplan-Meier estimator under independent censoring; the TMLE is
guaranteed to be within the bounds of the parameter space for θ. The latter property guarantees
that the estimated RMST for each study arm is nonnegative, which is not guaranteed for the AIPW
estimator; a negative RMST would be non-interpretable.
Using the R code in the Supplementary Materials, computation of the proposed estimator took
3.5 and 16.7 minutes for sample sizes 500 and 2000, respectively, for one of the simulated trials.
This computation time includes fitting the initial estimators for the event hazard, as well as the
censoring and the treatment mechanism. It was performed using R version 3.2.3 on a MacBook
Air with an Intel Core i5 1.3 GHz processor and 4GB of RAM.

19

Table 1: Simulation results for studies of size n = 500. The bias, variance (VAR), and mean
squared error (MSE) are displayed for the Kaplan-Meier (KM), unadjusted inverse probability
weighted (Unadj.IPW), adjusted IPW (IPW), augmented IPW (AIPW) and proposed targeted minimum loss based (TMLE) estimator. The relative MSE (RMSE) is the ratio of the MSE for the KM
estimator to the other estimators.
Non-informative censoring
Scenario

Estimator
KM
Unadj.IPW
A
IPW
AIPW
TMLE
KM
Unadj.IPW
B
IPW
AIPW
TMLE
KM
Unadj.IPW
IPW
C
AIPW
TMLE
Informative censoring
Scenario

A

B

C

Estimator
KM
Unadj.IPW
IPW
AIPW
TMLE
KM
Unadj.IPW
IPW
AIPW
TMLE
KM
Unadj.IPW
IPW
AIPW
TMLE

Zero Treatment Effect
Bias
Var
MSE RMSE
0.031 35.23 35.23 1.000
0.032 35.59 35.59 0.990
0.001 31.76 31.76 1.109
-0.020 31.22 31.21 1.128
-0.019 30.96 30.96 1.138
0.062 36.09 36.09 1.000
0.063 36.46 36.46 0.990
0.069 35.89 35.89 1.005
0.069 35.56 35.56 1.015
0.069 35.27 35.27 1.023
-0.048 35.52 35.52 1.000
-0.048 35.90 35.90 0.989
-0.046 36.41 36.41 0.976
-0.046 36.05 36.05 0.985
-0.046 35.76 35.76 0.993

Positive Treatment Effect
Bias
Var
MSE RMSE
0.013 26.54 26.54 1.000
0.044 26.86 26.86 0.988
0.004 24.29 24.29 1.093
-0.027 23.95 23.95 1.108
-0.074 23.74 23.74 1.118
-0.059 27.40 27.40 1.000
-0.029 27.72 27.72 0.988
-0.046 27.26 27.26 1.005
-0.087 26.94 26.94 1.017
-0.133 26.71 26.73 1.025
-0.035 26.91 26.91 1.000
-0.004 27.22 27.22 0.989
-0.013 27.50 27.50 0.978
-0.048 27.23 27.23 0.988
-0.092 27.00 27.01 0.996

Zero Treatment Effect
Bias
Var
MSE RMSE
0.802 31.48 32.12 1.000
0.794 31.61 32.25 0.996
-0.108 29.48 29.49 1.089
-0.045 29.22 29.22 1.099
-0.044 29.08 29.08 1.105
0.368 33.25 33.38 1.000
0.357 33.34 33.52 0.996
-0.084 33.68 33.68 0.991
-0.050 33.36 33.36 1.000
-0.050 33.09 33.09 1.009
0.039 31.88 31.88 1.000
0.025 32.02 32.02 0.996
0.022 32.66 32.66 0.976
0.041 32.51 32.51 0.980
0.041 32.24 32.24 0.989

Positive Treatment Effect
Bias
Var
MSE RMSE
1.760 22.48 25.58 1.000
1.775 22.59 25.73 0.994
0.079 22.92 22.93 1.116
0.097 22.44 22.44 1.140
0.045 22.18 22.18 1.153
0.882 23.12 23.90 1.000
0.893 23.23 24.03 0.994
0.046 24.31 24.31 0.983
0.056 24.13 24.13 0.990
0.004 23.93 23.92 0.999
0.039 24.02 24.02 1.000
0.046 24.13 24.13 0.995
0.053 24.90 24.90 0.965
0.051 24.77 24.77 0.970
0.001 24.65 24.65 0.974

20

Table 2: Simulation results for studies of size n = 2000. The bias, variance (VAR), and mean
squared error (MSE) are displayed for the Kaplan-Meier (KM), unadjusted inverse probability
weighted (Unadj.IPW), adjusted IPW (IPW), augmented IPW (AIPW) and proposed targeted minimum loss based (TMLE) estimator. The relative MSE (RMSE) is the ratio of the MSE for the KM
estimator to the other estimators.
Non-informative censoring
Scenario

Estimator
KM
Unadj.IPW
A
IPW
AIPW
TMLE
KM
Unadj.IPW
B
IPW
AIPW
TMLE
KM
Unadj.IPW
IPW
C
AIPW
TMLE
Informative censoring
Scenario

A

B

C

Estimator
KM
Unadj.IPW
IPW
AIPW
TMLE
KM
Unadj.IPW
IPW
AIPW
TMLE
KM
Unadj.IPW
IPW
AIPW
TMLE

Zero Treatment Effect
Bias
Var MSE RMSE
-0.069 8.86 8.86 1.000
-0.070 8.95 8.96 0.990
-0.057 8.04 8.04 1.102
-0.052 7.88 7.88 1.125
-0.051 7.81 7.81 1.134
0.022 8.92 8.92 1.000
0.022 9.01 9.01 0.990
0.019 8.77 8.77 1.018
0.017 8.69 8.69 1.027
0.017 8.62 8.62 1.035
-0.008 8.78 8.78 1.000
-0.008 8.87 8.87 0.990
-0.008 8.91 8.91 0.986
-0.009 8.82 8.82 0.996
-0.009 8.82 8.82 0.996

Positive Treatment Effect
Bias
Var MSE RMSE
-0.015 6.79 6.79 1.000
0.016 6.87 6.87 0.988
0.027 6.12 6.12 1.110
-0.015 6.00 6.00 1.133
-0.060 5.95 5.95 1.142
-0.060 6.68 6.68 1.000
-0.029 6.76 6.76 0.989
-0.033 6.60 6.60 1.012
-0.075 6.50 6.51 1.027
-0.119 6.45 6.46 1.035
-0.035 6.77 6.77 1.000
-0.004 6.85 6.85 0.989
-0.005 6.87 6.87 0.986
-0.048 6.80 6.80 0.996
-0.038 6.79 6.79 0.997

Zero Treatment Effect
Bias
Var MSE RMSE
0.797 7.78 8.42 1.000
0.790 7.82 8.44 0.997
-0.082 7.32 7.32 1.150
-0.022 7.20 7.20 1.168
-0.023 7.15 7.15 1.178
0.407 8.08 8.24 1.000
0.396 8.12 8.27 0.997
-0.054 8.15 8.15 1.011
-0.016 8.08 8.08 1.020
-0.016 8.02 8.02 1.028
-0.012 8.12 8.12 1.000
-0.025 8.16 8.16 0.995
-0.028 8.23 8.23 0.987
-0.009 8.19 8.19 0.992
-0.009 8.12 8.12 1.000

Positive Treatment Effect
Bias
Var MSE RMSE
1.646 5.65 8.36 1.000
1.661 5.68 8.44 0.991
-0.041 5.66 5.67 1.478
-0.012 5.53 5.53 1.513
-0.057 5.48 5.48 1.526
0.863 5.76 6.50 1.000
0.875 5.78 6.55 0.993
0.012 6.05 6.05 1.075
0.020 5.98 5.98 1.087
-0.024 5.93 5.93 1.097
-0.006 5.85 5.85 1.000
0.002 5.88 5.88 0.995
-0.010 5.98 5.98 0.978
-0.015 5.96 5.95 0.983
-0.061 5.90 5.91 0.991

21

8

Discussion

Under random censoring, our estimator is consistent if either the outcome or censoring model is
correctly specified, i.e., our estimator is doubly robust. This is in contrast to the proportional hazard
model, which relies exclusively on assumptions on the outcome model. When the dimension of
the baseline variables is large relative to the sample size, the curse of dimensionality precludes
the use of nonparametric estimators or saturated parametric models (Robins and Ritov, 1997). A
potential way to address this is to incorporate data-adaptive model selection in constructing the
initial estimators in step 1 of the TMLE procedure, such as model stacking (Wolpert, 1992) or
super learning (van der Laan et al., 2007). The asymptotic properties of the resulting estimator
then require conditions (f)-(g) in Theorem 3 of the Supplementary Materials. These conditions
would hold automatically for the MLE in a parametric model, but need to be verified for dataadaptive estimators. van der Laan (2014) proposed an estimator for the case of a cross-sectional
study that relaxes assumption (f), and we conjecture that this approach might be generalizable to
our problem.
In our presentation of the TMLE algorithm we have assumed that the initial estimators ĝR and
ĝA contain saturated terms for treatment, time, and their interaction. When using data-adaptive
methods, such restriction can be avoided by including the aforementioned saturated terms in the
logistic regression models (20)-(21) in step 2 of the TMLE algorithm.
Precision gains from adjustment for prognostic baseline variables can be converted into shorter
duration trials by using information monitoring. That is, the trial continues until a prespecified information level is achieved. Since improved precision implies (asymptotically) a faster information
accrual rate, the trial duration may be shorter.
Our method assumes that censoring is confounded with the time to event only by baseline
variables. In the presence of time dependent confounding between censoring and the event time,
our proposal may be adapted by augmenting the censoring and outcome models to include timevarying confounders. It may be possible to retain the properties (i)-(iv) in this context by extending
the proof techniques from the Supplementary Materials. These properties are still relevant under
time dependent confounding, since the use of proportional hazard models often yields biased estimators in this context, as discussed by Cole et al. (2003).
We defined independent censoring as C⊥⊥(Ta , W ) | A. Another possible definition of independent censoring, which is more commonly used when discussing the unadjusted estimator, is
C⊥⊥Ta | A (which leaves out W altogether). The latter assumption is weaker (less restrictive)
than the former. The covariate adjustment method of Zhang (2014) guarantees enhanced efficiency properties under the latter definition. Our TMLE estimator requires the former definition
to achieve the enhanced efficiency property (ii). On the other hand, our TMLE estimator has the
advantage over Zhang (2014) and θ̂km of remaining consistent under violations to the assumption
that C⊥⊥Ta | A, as long as random censoring holds and at least one of the censoring or survival
distributions is consistently estimated.
Most clinical research studies use discrete time scales to measure the time to event. If time
is measured on a continuous scale, implementation of our methods requires discretization. The
specific choice of the discretization intervals may be guided by what is clinically relevant. For
example, in clinical applications with time to death outcomes, the clinically relevant scale would
typically be a day. In the absence of clinical criteria to guide the choice of discretization level, a
concern is that too coarse of a discretization may lead to potentially meaningful information losses.
22

A question for future research is how to optimally set the level of discretization in order to trade
off information loss versus estimator precision. Another area for future research is to consider
discretization levels that get finer with sample size.

9

Acknowledgments

This work was supported by the Patient-Centered Outcomes Research Institute [ME-1306-03198],
the U.S. Food and Drug Administration [HHSF223201400113C], and the National Institute of
Neurological Disorders and Stroke (NINDS) [U01NS062851]. We thank the CLEAR III trial
team. This work is solely the responsibility of the authors and does not represent the views of the
above people and agencies.

A

Supplementary Materials
• Proofs for the main results in the paper. (proofs.pdf file)
• R functions to compute the proposed estimator. (functions.R file)

References
Heejung Bang and James M Robins. Doubly robust estimation in missing data and causal inference
models. Biometrics, 61(4):962–973, 2005.
P.J. Bickel, C.A.J. Klaassen, Y. Ritov, and J. Wellner. Efficient and Adaptive Estimation for Semiparametric Models. Springer-Verlag, 1997.
Weihua Cao, Anastasios A Tsiatis, and Marie Davidian. Improving efficiency and robustness of
the doubly robust estimator for a population mean with incomplete data. Biometrika, 96(3):
723–734, 2009.
Pei-Yun Chen and Anastasios A Tsiatis. Causal inference on the difference of the restricted mean
lifetime between two groups. Biometrics, 57(4):1030–1038, 2001.
Stephen R Cole and Miguel A Hernán. Adjusted survival curves with inverse probability weights.
Computer Methods and Programs in Biomedicine, 75(1):45–49, 2004.
Stephen R Cole, Miguel A Hernán, James M Robins, Kathryn Anastos, Joan Chmiel, Roger Detels,
Carolyn Ervin, Joseph Feldman, Ruth Greenblatt, Lawrence Kingsley, et al. Effect of highly
active antiretroviral therapy on time to acquired immunodeficiency syndrome or death using
marginal structural models. American Journal of Epidemiology, 158(7):687–694, 2003.
D. R. Cox. Regression models and life-tables (with discussion). Journal of the Royal Statistical
Society. Series B, 34(2):187–220, 1972.
Iván Dı́az, Elizabeth Colantuoni, and Michael Rosenblum. Enhanced precision in the analysis of
randomized trials with ordinal outcomes. Biometrics, 2015. ISSN 1541-0420.
23

B Efron et al. Bootstrap methods: Another look at the jackknife. The Annals of Statistics, 7(1):
1–26, 1979.
Susan Gruber and Mark J. van der Laan. Targeted minimum loss based estimator that outperforms
a given estimator. The International Journal of Biostatistics, 8(1):1–22, 2012.
Jinyong Hahn. On the role of the propensity score in efficient semiparametric estimation of average
treatment effects. Econometrica, pages 315–331, 1998.
Alan E Hubbard, Mark J Van Der Laan, and James M Robins. Nonparametric locally efficient
estimation of the treatment specific survival distribution with right censored data and covariates
in observational studies. In Statistical Models in Epidemiology, the Environment, and Clinical
Trials, pages 135–177. Springer, 2000.
Xiaomin Lu and Anastasios A Tsiatis. Semiparametric estimation of treatment effect with timelagged response in the presence of informative censoring. Lifetime Data Analysis, 17(4):566–
593, 2011.
Kelly L. Moore and Mark J. van der Laan. Increasing power in randomized trials with right
censored outcomes through covariate adjustment. Journal of Biopharmaceutical Statistics, 19
(6):1099–1131, 2009a.
Kelly L Moore and Mark J van der Laan. Covariate adjustment in randomized trials with binary outcomes: Targeted maximum likelihood estimation. Statistics in Medicine, 28(1):39–64,
2009b.
Layla Parast, Lu Tian, and Tianxi Cai. Landmark estimation of survival and treatment effect in a
randomized clinical trial. Journal of the American Statistical Association, 109(505):384–394,
2014.
J Pfanzagl and W Wefelmeyer. Contributions to a general asymptotic statistical theory. Statistics
& Risk Modeling, 3(3-4):379–388, 1985.
James M Robins and Ya’acov Ritov. Toward a curse of dimensionality appropriate (coda) asymptotic theory for semi-parametric models. Statistics in Medicine, 16(3):285–319, 1997.
J.M. Robins and A. Rotnitzky. Recovery of information and adjustment for dependent censoring
using surrogate markers. In AIDS Epidemiology, Methodological issues. Bikhäuser, 1992.
J.M. Robins, A. Rotnitzky, and L.P. Zhao. Estimation of regression coefficients when some regressors are not always observed. Journal of the American Statistical Association, 89(427):846–866,
September 1994.
Andrea Rotnitzky and James M Robins. Inverse probability weighting in survival analysis. Encyclopedia of Biostatistics, 2005.
Andrea Rotnitzky, Quanhong Lei, Mariela Sued, and James M Robins. Improved double-robust
estimation in missing data and causal inference models. Biometrika, 99(2):439–456, 2012.

24

Patrick Royston and Mahesh KB Parmar. The use of restricted mean survival time to estimate
the treatment effect in randomized clinical trials when the proportional hazards assumption is in
doubt. Statistics in Medicine, 30(19):2409–2421, 2011.
Donald B Rubin. Multiple Imputation for Nonresponse in Surveys. John Wiley & Sons, 1987.
Daniel O. Scharfstein, Andrea Rotnitzky, and James M. Robins. Adjusting for nonignorable dropout using semiparametric nonresponse models: Rejoinder. Journal of the American Statistical
Association, 94(448):pp. 1135–1146, 1999. ISSN 01621459.
Michael Schemper. Cox analysis of survival data with non-proportional hazard functions. Journal
of the Royal Statistical Society. Series D. (The Statistician), pages 455–465, 1992.
Ori M Stitelman, Victor De Gruttola, and Mark J van der Laan. A general implementation of tmle
for longitudinal data applied to causal inference in survival analysis. The International Journal
of Biostatistics, 8(1), 2011.
Zhiqiang Tan. A distributional approach for causal inference using propensity scores. Journal of
the American Statistical Association, 101(476):1619–1637, 2006.
Zhiqiang Tan. Bounded, efficient and doubly robust estimation with inverse weighting. Biometrika,
97(3):661–682, 2010.
Lu Tian, Lihui Zhao, and L. J. Wei. Predicting the restricted mean event time with the subject’s baseline covariates in survival analysis. Biostatistics, 15(2):222–233, 2014. doi:
10.1093/biostatistics/kxt050.
Anastasios A Tsiatis, Marie Davidian, Min Zhang, and Xiaomin Lu. Covariate adjustment for twosample treatment comparisons in randomized clinical trials: A principled yet flexible approach.
Statistics in Medicine, 27(23):4658–4677, 2008.
Mark J van der Laan. Targeted estimation of nuisance parameters to obtain valid statistical inference. The International Journal of Biostatistics, 10(1):29–57, 2014.
M.J. van der Laan and J.M. Robins. Unified Methods for Censored Longitudinal Data and Causality. Springer, New York, 2003.
M.J. van der Laan and D. Rubin. Targeted maximum likelihood learning. The International Journal
of Biostatistics, 2(1):Article 11, 2006.
M.J. van der Laan, E. Polley, and A. Hubbard. Super learner. Statistical Applications in Genetics
& Molecular Biology, 6(25):Article 25, 2007.
A. W. van der Vaart. Asymptotic Statistics. Cambridge University Press, 1998.
Jon A Wellner and Yihui Zhan. Bootstrapping z-estimators. University of Washington Department
of Statistics Technical Report, 308, 1996.

25

Elizabeth J Williamson, Andrew Forbes, and Ian R White. Variance reduction in randomised trials
by inverse probability weighting using the propensity score. Statistics in Medicine, 33(5):721–
737, 2014.
David H Wolpert. Stacked generalization. Neural Networks, 5(2):241–259, 1992.
Jun Xie and Chaofeng Liu. Adjusted kaplan–meier estimator and log-rank test with inverse probability of treatment weighting for survival data. Statistics in Medicine, 24(20):3089–3110, 2005.
Min Zhang. Robust methods to improve efficiency and reduce bias in estimating survival curves
in randomized clinical trials. Lifetime Data Analysis, 21(1):119–137, 2014. doi: 10.1007/
s10985-014-9291-y.
Min Zhang, Anastasios A Tsiatis, and Marie Davidian. Improving efficiency of inferences in
randomized clinical trials using auxiliary covariates. Biometrics, 64(3):707–715, 2008.
Lihui Zhao, Lu Tian, Hajime Uno, Scott D Solomon, Marc A Pfeffer, Jerald S Schindler, and
Lee Jen Wei. Utilizing the integrated difference of two survival functions to quantify the treatment contrast for designing, monitoring, and analyzing a comparative clinical study. Clinical
Trials, 9(5):570–577, 2012.
Lihui Zhao, Brian Claggett, Lu Tian, Hajime Uno, Marc A. Pfeffer, Scott D. Solomon, Lorenzo
Trippa, and L. J. Wei. On the restricted mean survival time curve in survival analysis. Biometrics,
72(1):215–221, 2016. ISSN 1541-0420. doi: 10.1111/biom.12384. URL http://dx.doi.org/10.
1111/biom.12384.

26

