Survival models and health sequences
Peter McCullagh∗July 3, 2014

Abstract
Survival studies often generate not only a survival time for each patient but
also a sequence of health measurements at annual or semi-annual check-ups
while the patient remains alive. Such a sequence of random length accompanied by a survival time is called a survival process. Ordinarily robust health isassociated with longer survival, so the two parts of a survival process cannot
be assumed independent. This paper is concerned with a general technique—reverse alignment—for constructing statistical models for survival processes.
A revival model is a regression model in the sense that it incorporates covariate and treatment effects into both the distribution of survival times and thejoint distribution of health outcomes. The revival model also determines aconditional survival distribution given the observed history, which describes
how the subsequent survival distribution is determined by the observed progression of health outcomes.

Keywords: interference; preferential sampling; quality-of-life; revival process; semirevival time; reverse alignment; stale values;

1

Survival studies

A survival study in one in which patients are recruited according to well-defined
selection criteria and their health status monitored on a regular or intermittent
schedule until the terminal event, here assumed to be fatal. Covariates such assex and age are recorded at the time of recruitment, and, if there is more than
one treatment level, the assignment is presumed to be randomized. In a simple
survival study, the health status Y (t) at time t is a bare-bones binary variable,
dead or alive, and the entire process is then summarized by the time T > 0 spentin state 1, i.e. the survival time. In a survival study with health monitoring, Y (t) isa more detailed description of the state of health or quality of life of the individual,
containing whatever information—pulse rate, cholesterol level, cognitive score orCD4 cell count—is deemed relevant to the study. The goal may be to study theeffect of treatment on survival time, or to study its effect on quality of life, or to
predict the subsequent survival time of patients given their current health history.
∗ Department of Statistics, University of Chicago, 5734 University Ave, Chicago, Il 60637, U.S.A.
E-mail: pmcc@galton.uchicago.edu

1Survival studies with intermittent health monitoring are moderately common,and likely to become more so as health records become available electronically for
research purposes. Within the past few years, several issues of the journal Lifetime
Data Analysis have been devoted to problems connected with studies of exactly
this type. For a good introduction, with examples and a discussion of scientificobjectives, see Diggle, Sousa and Chetwynd (2006), Kurland, Johnson, Eglestonand Diehr (2009) or Farewell and Henderson (2010). Section 8 of van Houwelingen
and Putter (2012) is recommended reading.
In practice, the patient’s health status is measured at recruitment (t = 0), andregularly or intermittently thereafter while the patient remains alive. To emphasize
the distinction between the observation times and observation values, each time
is called an appointment date, the set of dates is called the appointment schedule. Apart from covariate and treatment values, a complete uncensored observation
on one patient (T, t, Y [t]) consists of the survival time T > 0, the appointment
schedule t ⊂ [0, T ), and the health status measurements Y [t] at these times. To
accommodate patients whose record is incomplete, a censoring indicator variable isalso included. In that case, the censoring time is usually, but not necessarily, equal
to the date of the most recent appointment.
In the sense that the health status is measured over time on each patient, asurvival study is a particular sort of longitudinal study. Certainly, temporal andother correlations are expected and must be accommodated. But the distinguishingfeature, that each sequence is terminated by failure or censoring, gives survivalprocess models a very distinct character: as an absorbing state, death contradicts
stationarity. For a good survey of the goals of such studies and the modelingstrategies employed, see Kurland, Johnson, Egleston and Diehr (2009).
The goal of this paper is not so much to recommend a particular statisticalmodel, as to explore a general mathematical framework for the construction ofsurvival-process models, permitting easy computation of the likelihood function
and parameter estimates, and straightforward derivation of predictive distributions
for individual survival times. For example, the paper has nothing to say on the
choice between proportional hazards and accelerated lifetimes for accommodating
treatment effects. Apart from reservations concerning the use of time-evolving covariates, all standard survival models are acceptable within the framework. Nor has
the paper anything to contribute to the choice between Bayesian and non-Bayesian
methods of analysis; prior distributions are not discussed, so either approach can beused. Administrative complications of the sort that are inevitable in medical andepidemiological research will be ignored for the most part, so no attempt is madeto provide a complete turnkey package. For example, the paper has little to say
about how best to handle incomplete records other than to recognize that censoring and delayed reporting are issues that must be addressed—again using standard
well-developed methods. Since most of the computations needed for model fitting
and parameter estimation are relatively standard and need not involve specialized
Markov chain or Monte Carlo algorithms, detailed discussion of computational techniques is omitted. The emphasis is on statistical ideas and principles, strategies formodel formulation, sampling schemes, and the distinction between time-dependent
variables and time-evolving variables in the definition of treatment effects.
2

2

Joint modelling

The medical and biostatistical literature contains numerous examples of studies involving both successive measurements on each patient, such as CD4 lymphocyte
cell counts, together with survival time (Lagakos, 1976; DeGruttola and Tu, 1994;
Faucett and Thomas 1996; Guo and Carlin, 2004; Fieuws, Verbeke, Maes and Vanrenterghem, 2008). Geriatric studies seldom focus exclusively on survival time, but
tend to emphasize variables related to quality of life, such as overall physical andmental health, mobility, independence, memory loss, mental acuity, and so on. In
the statistical literature, survival studies with health monitoring are called longitudinal studies with time-to-event data (Wulfsohn and Tsiatis, 1997; Henderson, Diggle
and Dobson, 2000; Xu and Zeger, 2001; Tsiatis and Davidian, 2004; Rizopoulos,2012). Although there are variations in model formulation and implementation, allauthors are agreed on the need for a joint distribution covering both survival time
and the progression of health outcomes.
The joint modelling approach begins with a pair (η, T ) consisting of a latentrecurrent temporal process η(t) together with a positive random variable T < ∞.
This joint distribution determines the distribution of the observable processη(t) t < T
Y (t) =[t≥T
by restriction to (−∞, T ), or censoring at T . The process constructed in this way
has the same domain t ∈ < as the unobservable recurrent process, but the state
space includes an additional absorbing value, here labelled [. Equivalently, butslightly less conveniently, Y may be defined as a process on the random domain(−∞, T ) or (0, T ). In either case the survival time T is a function of Y .
In practice, the recurrent process usually has a fairly simple form such as astationary Gaussian process, and T is determined by a Cox model whose intensity
is a latent process not independent of η. The values (ηi , Ti ) are assumed independentfor distinct individuals, and identically distributed for individuals sharing the same
covariate values. Regardless of how the survival process is constructed, the joint
distribution is determined by its finite-dimensional restrictions. In other words, weneed to specify for each finite subset s ⊂ [0, ∞), the joint distribution of Y [s], orthe joint distribution of (T, Y [s]), or, equivalently, the joint distribution of (T, Y [t]),where t = s ∩ [0, T ). For a more detailed account of specific processes, see chapters
13–16 of Fitzmaurice, Davidian, Verbeke and Molenberghs (2009).
In order to avoid some of the technical difficulties associated with joint modelling the suggestion put forward in this paper is to approach the problem froma different angle—literally in reverse. Reverse alignment avoids the intermediate
recurrent process by considering probability distributions for the non-recurrent observable process Y directly. It is discussed as one of several options in Table 2of Kurland, Johnson, Egleston and Diehr (2009), and is mentioned in section 8.3
of van Houwelingen and Putter (2012), so the idea is not new. This paper explores the probabilistic and statistical implications of time-reversal, with a focus
on exchangeability and distributional factorization for the entire survival process.

3

The implications for sampling, the consequences for survival prediction, and theinterpretation and estimation of treatment effects are also considered.

3

Reverse alignment

3.1

The survival process

A survival process Y is a stochastic process defined for real t, in which Yi (t) is
the state of health or quality of life of patient i at time t, usually measured fromrecruitment. In a simple survival process, the state space R = {0, 1} is sufficient toencode only the most basic of vital signs, dead or alive; more generally, the state
space is any set large enough to encode the observable state of health or quality of lifeof the patient at one instant in time. Flatlining is the distinguishing characteristicof a survival process, i.e. [ ∈ R is an absorbing state such that Y (t) = [ implies
Y (t0 ) = [ for all t0 ≥ t. The survival time is the time to failure:
Ti = sup{t : Yi (t) 6= [};t≥0

it is presumed that Yi (0) 6= [ at recruitment, so Ti > 0. This definition is quite
general, and does not exclude immortality, i.e. T = ∞ with positive probability. In
all of the models considered here, however, survival time is finite with probability
one.

3.2

Administrative and other schedules

Since the appointment schedule is a random subset t ⊂ [0, T ), it is obviously informative for survival: T > max(t). If better health is associated with longer survival,
we should expect patients who are initially frail to have shorter health records
than patients who are initially healthy. In other words, even if the trajectories for
distinct individuals may be identically distributed, the first component of a short
health-status record should not be expected to have the same distribution as thefirst component of a longer record. On the contrary, any model such that record
length is independent of record values must be regarded as highly dubious for survival studies. It is necessary, therefore, to address the nature of the information
contained in t.
Consider a patient who has had appointments on k occasions t(k) = (t0 < · · · <tk−1 ). The sequence Y [t(k) ] of recorded health values may affect the scheduled
date tk for the next appointment: for example, patients in poor health needing
more careful monitoring may have short inter-appointment intervals. Whatever thescheduled date may be, the appointment is null unless tk < T . The assumption
used in this paper is sequential conditional independence, namely thattk ⊥⊥ Y | (T, t(k) , Y [t(k) ]).

(1)

In other words, the conditional distribution of the random interval tk − tk−1 may
depend on the observed history Y [t(k) ], but not on the subsequent health trajectory4except through T . Here, tk may be infinite (or null) with positive probability, inwhich case the recorded sequence is terminated at tk−1 .
The schedule is said to be administrative if tk is a deterministic function of
the pair (t(k) , Y [t(k) ]), implying that the conditional distribution (1) is degenerate.
Eventually, for some finite k, the patient dies or is censored at time T ∈ (tk−1 , tk )while the next appointment is pending, so the recorded schedule is t = t(k) =t(k+1) ∩ [0, T ). Equivalently, the last recorded value is (tk , [).
While the sequential conditional independence assumption is mathematically
clear-cut, the situation in practice may be considerably more muddy. Consider,for example, the CSL1 trial organized by the Copenhagen Study Group for Liver
Diseases in the 1960s to study the effect of prednizone on the survival of patientsdiagnosed with liver cirrhosis. In this instance Y (·) is a composite blood coagulation index called the prothrombin level: details can be found in Andersen, Hansen
and Keiding (1991). Beginning at death, the reverse-time mean intervals between
appointments are 77, 210 and 251 days, while the medians are 21, 165 and 292 days.
In other words, half of the patients who died had their final appointment within
the last three weeks of life. It is evident that the appointment intensity increases
as s → 0 in reverse time, which is not, in itself, a violation of (1). However, one
might surmise that the increased intensity is related to the patient’s state of healthor perception thereof. Condition (1) implies that the appointment intensity does
not depend on the blood coagulation index other than at earlier appointments, and
it is then unclear to what extent the condition may be violated by patient-initiatedappointments. Liestøl and Andersen (2002, section 4.1) note that 71 off-schedule
appointments occurred less than 10 days prior to death, the majority of which were
patient-initiated. They also examine the effect on hazard estimates of excluding
unscheduled prothrombin measurements.
Although we refer to Y (·) generically as the patient’s state of health, this description is not to be taken literally. The actual meaning depends on what has infact been measured: in general, Y (·) is only one component or one aspect of patient
health.

3.3

The revival process

On the assumption that the survival time is finite, the time-reversed process
Zi (s) = Yi (Ti − s)is called the revival process. Thus, Zi (s) is the state of health of patient i at times prior to failure, and Zi (Ti ) = Yi (0) is the value at recruitment. By construction,Z(s) = [ for s < 0, and Z(s) 6= [ for s > 0. Although Z is defined in reverse time,
the temporal evolution via the survival process occurs in real time: by definition,
Z(·) is not observable component-wise until the patient dies.
The transformation Y 7→ (T, Z) is clearly invertible; it may appear trivial, and ina sense it is trivial. Its one key property is that the revival process Z and the randomvariable T are variation independent. In the statistical models considered here,
variation independence may be exploited through the revival assumption, which

5states that the revival process and the survival time are statistically independent.
More generally, Z ⊥⊥ T | X if covariates are present. The revival assumption
provides a convenient starting point for model construction; it is not critical to any
part of the theory presented here.
Table 1: Average prothrombin levels indexed by T and t.
Survival
Time t after recruitment (yrs)time (T ) 0–1 1–2 2–3 3–4 4–5 5–6 6–7 7–8
0–1
1–2
2–3
3–4
4–5
5–6
6–7
7–8
8+

58.0
72.5
72.6
69.8
68.5
70.5
81.8
84.4
77.3

66.4
73.2
71.2
75.7
77.3
73.6
88.8
73.6

66.0
68.5
72.5
73.5
81.1
88.1
87.0

54.2
74.6
57.1
80.6
92.1
74.1

57.7
64.5
79.4
85.2
92.0

60.9
75.5
81.2
80.3

75.8
84.3
89.2

88.1
79.4

8+

84.7

The chief motivation for time reversal has to do with the effective alignment ofpatient records for comparison and signal extraction. Are the temporal patterns
likely to be more similar in two records aligned either by patient age or by recruitment date, or are they likely to be more similar in records aligned by reverse age(time remaining to failure)? Ultimately, the answer must depend on the context,
but the context of survival studies suggests that the latter may be more effective
than the former. Table 1 shows the averaged Y -values indexed by T and t for theprothrombin example discussed in more detail in section 6. It should be borne inmind that each cell is the average of 8–266 non-independent high-variability measurements, the larger counts occurring in the upper left cells. Alignment by reversetime is equivalent to counting leftwards from the main diagonal. Despite certain
anomalies in the table of averages, e.g. row 6, column 4, it is clear that reverse-time
is a more effective way of organizing the data to display the main trends in themean response: the forward- and reverse-time sums of squares (equally weighted)
are 543.0 and 1132.8 respectively, both on eight degrees of freedom.
Further confirmation is provided in Table 2, which shows the output from astandard equally-weighted analysis of variance applied to the table of averages,
with three factors, row, column and diagonal (reverse time), denoted by R, C and D
respectively. Compared with the residual mean square of 23.7, there is considerable
excess variation associated with rows (117.0) and with the reverse-time factor (77.8),
but not so much with columns (34.1). In other words, the means in Table 1 areexpressible approximately as αT + βT −t . Figures. 8.3 and 8.4 of van Houwelingen
and Putter (2012), which are not substantially different from Fig. 2 of this paper,
offer strong confirmation of this viewpoint in one further survival study involving
white blood cell counts for patients suffering from chronic myeloid leukemia. Foran application unrelated to survival, see example B of Cox and Snell (1981).
Model construction by reverse alignment may appear peculiar and unnatural in6

Table 2: ANOVA decomposition for Table 1
U/V
kPU Y k2 − kPV Y k2
(R + C + D)/(R + C)
544.4
(R + C + D)/(R + D)
238.9
(R + C + D)/(C + D)
818.8
RC/(R + C + D)
498.3

d.f.
7
7
7
21

M.S.
77.8
34.1
117.0
23.7

biological work, where the accepted wisdom is that an effect such as death cannot
precede its supposed causes, such as ill health. However, the contrarian viewpoint(Liestøl and Andersen (2002, section 4.0)), that proximity to death is the chief
cause of ill health, seems neither less compelling nor more helpful. The author’s
attitude is that metaphysical discussion along such lines is seldom productive andbest avoided.

3.4

Exchangeability

In the presence of covariates such as sex, age at recruitment or treatment status,
exchangeability is understood in the sense of McCullagh (2008, section 2) i.e. it
applies to each subset of patients having the same covariate value. For any such setof patients, it implies that the survival times Ti1 , . . . , Tin are identically distributed,
the record lengths #ti1 , . . . , #tin are identically distributed, the health-status variables Yi1 (t), . . . , Yin (t) are identically distributed, and likewise for the revival values
at any fixed revival time s. In this paper, therefore, baseline health status is the firstcomponent of Y, not a covariate. This is unavoidable for revival models: it is immaterial that Y (0) is measured prior to randomization and treatment assignment.
Exchangeability does not imply that Yi (0) is independent of the record length #ti .
Assuming that the record is complete, the observation for one patient consistsof a survival time T , a finite appointment schedule t ⊂ [0, T ), and a sequence oflength #t taking values in the space of medical records, here denoted by R. For
simplicity of notation in what follows, it is assumed that the appointment date isincluded in R. Then the sample space for the observation (T, Y [t]) on one patient
is
∞
[
Rk
S = (0, ∞) ×
k=0

in which the second component is the space of finite-length R-valued sequences.
For n patients i1 , . . . , in , the observations are independent if the joint distribution on S n factors in the usual way:
Pi1 ,...,in (A1 × · · · × An ) =nY

Pij (Aj ),j=1

where each Pi is a probability distribution on S, and A1 , . . . , An ⊂ S are arbitrary
events. In that circumstance, it is sufficient to describe the marginal distributions
Pi on S, which may depend on covariates xi . The observations on patients are7

infinitely exchangeable if Pi1 ,...,in is the marginal distribution of Pi1 ,...,in ,in+1 , andall joint distributions are unaffected by permutation of patients, who are always
assumed to be distinct individuals.
The implications of exchangeability are the same whether the record for eachpatient is expressed in terms of the survival process or the revival process. It implies
that the revival processes are identically distributed. Together with the revival
assumption, that Z and T are independent, it implies that Zi1 (s), . . . , Zin (s) are
identically distributed independently of the survival times Ti1 , . . . , Tin .

3.5

Covariates

In the absence of specific information to the contrary, responses for distinct units
are presumed to be identically distributed. In the great majority of situations,specific information does exist in the form of covariates or classification variables
or relationships. A covariate is a function i 7→ xi on the units, in principle known
for all units whether they occur in the sample or not. A covariate implies a specificform of inhomogeneity such that equality of covariates implies equality of response
distributions: xi = xj implies Yi ∼ Yj . In practice, approximate equality of x-values
also implies approximate equality of distributions. Likewise, a relationship is afunction on pairs of units such that R(i, i0 ) = R(j, j 0 ) implies (Yi , Yi0 ) ∼ (Yj , Yj 0 ) for
distinct pairs i 6= i0 , j 6= j 0 , provided that the two pairs also have the same covariate
values: (xi , xi0 ) = (xj , xj 0 ). Geographic distance and genetic distance are two
examples of symmetric relationships. The overarching principle is that differencesin distribution, marginal or joint, must be associated with specific inhomogeneities
in the experimental material.
The status of certain variables in specific survival studies may appear genuinely
unclear. The conventional rationalization, in which certain variables used for prediction are notionally ‘fixed’ or non-random and treated as covariates, is not especiallyhelpful for survival studies. Consider, for example, marital status as one variable
in a geriatric study in which the goal is to study both quality of life and survival
time. However it is defined, quality of life is a multi-dimensional response, a combination of mobility, independence, optimism, happiness, family support, and soforth. Marital status is a temporal variable known to be associated with survival
and with quality of life; one goal may be to predict survival given marital status,
or even to recommend a change of status in an effort to improve the quality of life.
Another example of a similar type is air quality and its relation to the frequency
and severity of asthmatic attacks (Laird, 1996). Should such a variable be regarded
as a covariate or as one component of the response? For survival studies, andfor longitudinal studies generally, the answer is very clear and very simple: every
time-evolving variable is necessarily part of the response process.
By definition, a temporal variable x is a function defined for every t ≥ 0. A
temporal variable is a covariate if it is also a function on the units, meaning thatthe entire function is determined and recorded at baseline. Usually this meansthat x is constant in time, but there are exceptions such as patient age: see also
section 3.6. Marital status and air quality, however, are not only temporal variables,
but variables whose trajectories evolve over real time; neither is available as a8

covariate at baseline.
With marital status as a component of the survival process, the joint distribution
may be used to predict the survival time beyond t of an individual whose marital
history and other health-status measurements at certain times prior to t are given.
For that purpose, it is necessary to compute the conditional distribution of T, or
more generally of Y, given the observed history Ht at the finite set of appointments
prior to t. For such calculations to make mathematical sense, marital status must
be a random variable, a function of the process Y. Thus, the statement ‘marital
status is a random process’ is not to be construed as a sociological statement about
the fragility of marriage or the nature of human relations; it is merely a mathematical assertion to the effect that probabilistic prediction is not possible without
the requisite mathematical structure of σ-fields Ht ⊂ Ht0 for t ≤ t0 and probability
distributions.

3.6

Treatment

Treatment refers to a scheduled intervention or series of interventions in which, atcertain pre-specified times following recruitment, the prescription for patient i is
switched from one arm to another. Thus, ai (t) is the treatment arm scheduled forpatient i at time t ≥ 0. In general, but crucially for revival models, a null level isneeded for t ≤ 0, including the baseline t = 0. The entire temporal trajectory ai (t)for t > 0 is determined by randomization and recorded at baseline. It does not evolve
over real time in response to the doctor’s orders or the patient’s perceived needs, so
it is not a time-evolving variable. Ordinarily, the random variables a1 (·), . . . , an (·)
are not independent. In the sense that it is recorded at baseline, ai (·) is a covariate;in the sense that it is a temporal function, it is a time-dependent covariate.
Apart from crossover trials, the distribution of a(·) is such that a switch oftreatment arms occurs only once, and then only immediately after recruitment.
Nonetheless, more general formulation is retained to underline the fact that treatment is a scheduled intervention such that ai (t) 6= ai (0), and thus not constant intime. Unlike the survival process, the treatment schedule does not evolve randomly
in real time.
It should be understood that a treatment arm is a protocol, also called a dynamic
treatment regime (Murphy, 2003), specifying the drug type, dose level, frequency,
manner of ingestion, and even the next appointment date, as a function of the current medical circumstances and health history. A treatment arm is not necessarily
associated with a specific drug or a specific medical therapy. Consider a hypertension study in which blood pressure Y (·) in conventional units is measured at regularappointments. Thus, one blue pill to be taken three times daily while blood pressure
exceeds 180, and one white pill twice daily if the pressure is between 160 and 180
is a treatment arm in which the actual therapy and dose level at time t0 dependon the outcome Y (t) at the most recent appointment t ≤ t0 . Another treatment
arm might reverse the colours or adjust the doses or change the therapy or reducethe inter-appointment interval if preceding Y -values are not encouraging. In thissetting, ai (t) denotes the assigned treatment arm, not the current drug or the dose,
so ai (·) is ordinarily constant for t > 0, even for dynamic chemotherapy strategies
9such as those discussed by Rosthøj, Keiding and Schmiegelow, (2012). For a compliant patient, the drug type and ingested dose can, in principle, be determined fromthe treatment arm and previously recorded Y -values. The key point is that eachtreatment arm be fully specified in advance, and the assignment be randomized atrecruitment.
Let āi (s) = ai (Ti − s) be the treatment arm expressed in revival time, so that,in the standard setting, āi (s) is null for s ≥ Ti . It is automatic that that Z ⊥⊥T | ā, because T is a function of ā. In the case of treatment, however, the crucial
assumption is lack of interference, i.e. the treatment assigned to one individual
has no effect on the response distribution for other individuals, and the treatment
protocol at one point in time has no effect on the response distribution at othertimes. For the latter, the statement is as follows. For each finite subset s ⊂ <+ ,the conditional distribution of Z[s] given the treatment schedule and survival timedepends only on the treatment arms ā[s] prevailing at the scheduled times, i.e.,Z[s] ⊥⊥ ā | ā[s].
For crossover trials in particular, this is a strong assumption denying carry-overeffects from earlier treatments or later treatments. It implies in particular that
Z(s) ⊥⊥ T | ā(s), which is primarily a statement about the one-dimensional marginal
distributions. Note, however, that the interference assumption is relatively benign
if ai (t) is constant for t > 0, as is ordinarily the case.
It is common practice in epidemiological work for certain time-evolving variables
to be handled as covariates, as if the entire trajectory were recorded at baseline.
This approach is perfectly reasonable for an external variable such as air quality
in an asthma study where lack of cross-temporal interference might be defensible.
It has the advantage of leading to simple well-developed procedures for effect estimation using marginal moments (Zeger and Liang, 1986; Zeger, Liang and Albert,
1988; Laird, 1996; Diggle, Heagerty, Liang and Zeger, 2002). The same approach isless convincing for an evolving variable such as marital status in a survival study, because the entire trajectory—suitably coded for t > Ti —would often contain enoughinformation to determine the survival time.

4

Survival prediction

4.1

Conditional distribution

Consider the simplest model in which observations for distinct patients are independent and identically distributed. To simplify matters further, problems related toparameter estimation are set aside. In other words, the survival time is distributed
according to F , and the revival processes given T = t is distributed as G(· | t). Given
the joint distribution, we are free to compute whatever conditional or marginal distribution is needed to address the inferential target.
We consider here the question of how the partial trajectory of Y affects thesubsequent survival prognosis. However the question is phrased, any suggestionof a causal mechanism is unwarranted: ultimately the answer is determined by
10

the appropriate conditional distribution. The problem is to predict the survival
time of an individual given the survival process Y [t(k) ] at the first k appointmentst(k) = (t0 < · · · < tk−1 ).
For positive real numbers s = (s1 > · · · > sk ), let gk (z; s | t) be the conditional
joint density given T = t of the health-status values
Z[s] = (Z(s1 ), . . . , Z(sk )) = (Y (T − s1 ), . . . , Y (T − sk )).
Under the conditional independence assumption (1), which implies non-preferential
appointment dates in the sense of Diggle, Menezes and Su (2010), the joint densityof (T, t(k) , Y [t(k) ]) at (t, t(k) , y) is a product of three factors:Yf (t) ×p(tj , yj | Hj , T = t)j<k= f (t) ×Yp(yj | Hj , T = t) ×j<kYp(tj | Hj , T = t)j<k= f (t) × gk (y; t − t(k) | t) ×Yp(tj | Hj , T = t),

(2)

j<k

where f = F 0 is the survival density, and Hj is the observed history (t(j) , Y [t(j) ])at time tj−1 . Without further assumptions, all three factors depend on t, meaningthat all three components are informative for survival prediction.
In all subsequent discussion concerning prediction, it is assumed that the appointment schedule is uninformative for prediction in the sense thatp(tk | Hk , T = t) = p(tk | Hk , T = ∞)(3)

for tk−1 < tk < t. This means that the next appointment is scheduled as if T = ∞,but it is not recorded unless tk < T . With this assumption, the third factor in (2)is constant in t and can be ignored. In other words, the distribution of the timeto the next scheduled appointment may depend on the patient’s medical history,
but is independent of the patient’s subsequent survival. Ordinarily, the scheduled
appointment is included as a component of the patient’s record only if it occurs in[0, T ) while the patient lives, implying that the partial appointment schedule t(k)is uninformative for subsequent survival. In particular, an administrative schedule
is uninformative.
A simple numerical example illustrates the idea. Suppose T is exponentiallydistributed with mean 10 years, and the revival process for s > 0 is a real-valued
Gaussian process with mean E(Z(s)) = βs/(1 + s) and covariance function δss0 +exp(−|s − s0 |) for s, s0 > 0. The observed health-status values at t = (0, 1, 2, 3) arey = (6.0, 4.5, 5.4, 4.0).
For β = 0, the conditional density is such that T −3 is exponential with mean 10;
the conditional density is shown for various values 0 ≤ β ≤ 2 in the left panel ofFig. 1, and for 4 ≤ β ≤ 8 on the right. Evidently, the conditional distribution
depends on both the observed outcomes and on the model parameters: the median

11

0.5
0.02

0.1

0.04

0.2

0.06

0.3

0.08

0.4

0.10

beta=8
beta=0.0

beta=40.0

0.00

beta=2.0

5

10

15

20

25

5

10

15

20

25

Figure 1: Conditional density of survival time for various values of β.
residual lifetime is not monotone in β. In applications where β is estimated with appreciable uncertainty, the predictive distribution is a weighted convex combinationof the densities illustrated.
The conditional survival distribution given Y [t(k) ] depends not only on the current or most recent value, but linearly on the entire vector. In particular, theconditional distribution does not have the structure of a regression model in which
the longitudinal variable enters as a time-dependent covariate without temporal interference. Thus, on the assumption that the joint model is adequate, issues related
to covariate confounding do not arise.

4.2

A simple Gaussian revival process

Under assumptions (1) and (3), the ratio of the conditional survival density at t
to the marginal density is proportional to the factor gk (y; t−t(k) ), in which y, t(k) are
fixed, and t the variable. This modification factor—the Radon-Nikodym derivative—
depends only on the revival process, not on the distribution of survival times. On
a purely mathematical level, it is precisely the likelihood function in the statisticalmodel for the k-dimensional variable Y [t(k) ] whose conditional distribution given
T = t is Gk (y; t − t(k) | t) for some value of the temporal offset parameter t > tk−1 .
Although not realistic for most applications, suppose that G is Gaussian withmean µ(s) = α + βs independent of t and linear in reverse time, and covariancefunction cov(Z(s), Z(s0 ) | t) = K(|s − s0 |). Then the log density ratio factor
log gk (y; t − t(k) | t) = const − 12 (y − µ[t − t])0 K −1 (y − µ[t − t]),is quadratic in t for t > tk−1 . After substituting α + β(t − t) for the mean function,and expressing the log density ratio as a quadratic in t, it can be seen that thepredictive density ratio at t > tk−1 is the density at βt of the Gaussian distribution
with mean−α + 10 K −1 (y + βt(k) ) (10 K −1 1) = ȳ − α + β t̄and variance 1/(10 K −1 1), where K has components K(ti − tj ). Ignoring the dependence on the data that comes from parameter estimation, the dependence of
12the predictive density ratio on the data for one patient comes through the weighted
averagesȳ = 10 K −1 y/(10 K −1 1),t̄ = 10 K −1 t/(10 K −1 1)(4)
for this particular individual.

4.3

Exchangeable Gaussian revival process

In a more natural Gaussian model, the revival processes for distinct patients are exchangeable but not necessarily independent. Revival models have much in common
with growth-curve models (Lee, 1988, 1991) in which Zi (s) = µ(s)+η0 (s)+ηi (s) is asum of two independent zero-mean Gaussian processes, and the mean function µ(s)is constant across individuals. Usually the common deviation η0 (·) is moderately
smooth but not stationary, perhaps fractional Brownian motion with η0 (0) = 0.
The idiosyncratic deviations are independent and identically distributed and they
incorporate measurement error, so ηi (·) is ordinarily the sum of a continuous process
and white noise. Thus, the Gaussian process is defined byE(Zis ) = µ(s)(5)

cov(Zis , Zi0 s0 ) = K0 (s, s0 ) + δii0 K1 (s, s0 ) + σ 2 δii0 δss0
for some suitable covariance functions K0 , K1 , each of which can be expected to havea variance or volatility parameter and a range parameter. In the case of fractional
Brownian motion, for example, K(s, t) ∝ sν + tν − |s − t|ν for some 0 < ν < 2,which governs the degree of smoothness of the random function.
For a new patient such that Y [t(k) ] = y, the conditional survival density pr(T ∈dt | data) given the data, including the outcomes for the new patient, is computed in
the same way as above. The second factor in (2) is the density at the observed outcomes of the Gaussian joint distribution whose means and covariances are specified
above. This involves all n + 1 patients.

4.4

Illustration by simulation

Figure 2 shows simulated data for 200 patients whose survival times are independent
exponential with mean five years. While the patient lives, annual appointments are
kept with probability 5/(5 + t), so appointment schedules in the simulation are
not entirely regular. Health status is a real-valued Gaussian process with meanE(Z(s)) = 10 + 10s/(10 + s) in reverse time, and covariancescov(Z(s), Z(s0 )) = (1 + exp(−|s − s0 |/5) + δss0 )/2
for s, s0 > 0, so there is an additive patient-specific effect in addition to temporalcorrelation. Values for distinct patients are independent and identically distributed.
This distribution is such that health-status plots in reverse time aligned by failure
show a stronger temporal trend than plots drawn in the conventional way. The
state of health is determined more by time remaining before failure than time since
recruitment. These trends could be accentuated by connecting successive dots for
13

20
8

10

12

14

16

18

20
18
16
14
12
10
8
0

5

10

15

20

25

30

−30

−25

−20

−15

−10

−5

0

Figure 2: Simulated health status sequences aligned by recruitment time (left) andthe same sequences aligned by failure time (right)
each individual, as in Fig. 2 of Sweeting and Thompson (2011), but this has notbeen done in Fig. 2.
Since the survival times are exponential with mean five, independent of covariates and treatment, the root mean squared prediction error using covariates onlyis five years. For fixed k ≥ 2, and a patient having at least k appointments, theconditional survival distribution given the first k health-status values has a standard deviation depending on the observed configuration, but the average standard
deviation is about 2.5 years, and the root mean squared prediction error is about 2.7
years. For this setting, the longitudinal variable is a reasonably effective predictor ofsurvival, and the prediction error is almost independent of k in the range 2–5. Thissummary does not tell the full story because certain y-configurations lead to very
precise predictions whereas others lead to predictive distributions whose standard
deviation exceeds five years.
The parameter settings used in this simulation may not be entirely representativeof the range of behaviours of the conditional survival distribution given Y [t]. If theratio of the between-patient to within-patient variance components is increased, theaverage variance of the conditional survival distribution decreases noticeably withk. For such settings, prediction using the entire health history is more effective than
prediction using the most recent value.

4.5

Recurrent health-related events

In certain circumstances the health outcome Y is best regarded as a point process,recording the occurrences of a specific type of non-fatal event, such as epileptic
or asthmatic attacks or emergency-room visits. In other words, Yi ⊂ < is the setof times at which patient i experiences the event. Then t = (0, tk ) is a boundedinterval, and the observation Y [t] = Y ∩ t is the set of events that occur between recruitment and the most recent appointment. This observation records theactual date of each event, which is more informative than the counting process#Y [(0, t1 )], . . . , #Y [(0, tk )] evaluated at the appointment dates. If there are recurrent events of several types, Y is a marked point process, and Y [t] is the set of all
14

events of all types that occur in the given temporal interval. The paper Schaubel
and Zhang (2010) is one of several papers in the October 2010 issue of Lifetime
Data Analysis, which is devoted to studies of this type.
In this situation, the frequency of the recurrent event may be constant over time,or it may vary in a systematic way. For example, the frequency may increase slowly
but systematically as a function of either age or time since recruitment. Alternatively, the frequency may be unrelated to age at recruitment, but may increase in
the last year of life as death approaches. In the former case, alignment of records
by failure time is ineffective; in the latter case, the revival processes for different
individuals have a common pattern, and alignment by failure tine is an effective
device for exploiting this commonality.
We consider here only the simplest sort of recurrent-event process in which therevival process is Poisson, there is a single event type, and the subset Y ∩ t = yof observed event times is finite. The mean measure of the revival process is Λ,which is non-atomic with intensity λ on the positive real line. The density ratio att > sup(t) is the probability density at the observed event configuration t − y as asubset of the reverse-time interval t − t, i.e.,Yg(y; t − t) = exp(−Λ(t − t))λ(t − y).
y∈y

x

x

x

x

x

x

x

x

x

x

0.4

0.6

exp(l0)

0.8

1.0

1.2

1.4

In particular, if the intensity is constant for s > 0, the density ratio is constant,and the event times are uninformative for survival. In other words, it is the temporal variation of the intensity function that makes the observed configuration y
informative for patient survival.

0

5

10

15

t

Figure 3: Likelihood functions (solid lines) for three point configurations, withpredictive hazard ratios (dashed lines)
For a specific numerical example, let λ(s) = (2 + s2 )/(1 + s2 ) be the revival
intensity, and let t = (0, 2) be the observation window. The revival intensity,
monotone decreasing with an asymptote of one, implies that the recurrent events aremoderately common at all ages, but their frequency increases as failure approaches.
Figure 3 shows the likelihood as a function of t ≥ 2 for three event configurations,15

y0 = ∅, y1 = {0.5, 1.2} and y2 = {0.2, 1.3, 1.9}. Since the likelihood function
is defined only up to an arbitrary multiplicative constant, the curves have beenadjusted so that they are equal at t = 20, or effectively at t = ∞. In place ofthe predictive survival distributions, we show instead the ratio of the predictive
hazard functions to the marginal hazards as dashed lines on the assumption thatthe marginal failure distribution is exponential with mean 5. Because of the formof the revival intensity, which is essentially constant except near the origin, thepredictive hazard functions are very similar in shape to the likelihood functions.

5

Parameter estimation

5.1

Likelihood factorization

The joint density for the observations in a revival model factors into two parts, one
involving only survival times, the other involving only the revival process. Moregenerally, if the revival assumption fails, the second factor is the conditional distribution of the revival process given T = t, so both factors depend on t. Althoughboth factors may involve the same covariates and treatment indicators, the parameters in the two parts are assumed to be unrelated, i.e. variation independent.
Thus the likelihood also factors, the first factor involving only survival parameterssuch as hazard modifiers associated with treatment and covariates, the second factor involving only health-status parameters such as temporal trends and temporalcorrelations. In other words, the two factors can be considered separately and independently, either for maximum likelihood estimation or for Bayesian operations.
The first stage in parameter estimation is to estimate the survival distribution
F together with treatment and covariate effects if needed. Whether the model forsurvival times is finite-dimensional or infinite-dimensional, this step is particularly
simple because the first factor involves only the survival times and survival distribution. The standard assumption of independent survival times for distinct patients
simplifies the problem even further. Exponential, gamma and Weibull models are
all feasible, with treatment effects included in the standard way.
For the Cox proportional-hazards model, the situation is a little more complicated. First, the survival Rtime is finite with probability one if and only if the∞integrated hazard Λ(<+ ) = 0 λ(t) dt is infinite, which is not satisfied at all parameter points in the model. Second, the proportional-hazards likelihood function
depends only on baseline hazard values λ(t) in the range 0 ≤ t ≤ Tmax , where
Tmax is the maximum observed survival time, censored or uncensored. Thus, thelikelihood does not have a unique maximum, but every maximum has the property
that λ̂(t) = 0 for all 0 ≤ t ≤ Tmax except for failure times, at which λ̂ has a discreteatom. By common convention (Kaplan and Meier 1958; Cox 1972, §8) λ̂(t) = 0 for
t > Tmax , but this choice is not dictated by the likelihood function. Since the revival
model requires survival times to be finite with probability one, it is essential to restrict the space of hazards to those having an infinite integral, which rules out thestandard convention for λ̂. Equivariance under monotone temporal transformations
points to a mathematically natural choice λ̂(t) = ∞ for t > Tmax ; a less pessimistic

16option is to use a finite non-zero constant such as
λ̂(t) =

total number of failures
total person time at risk

for t > Tmax . Both of these maximize the proportional-hazards likelihood function—
restricted or unrestricted—and either one may be used in the revival model.
The second stage, which is to estimate the parameters in the revival process, isalso straightforward, but only if all records are complete with no censoring. Serial
dependence is inevitable in a temporal process, and there may also be independent
persistent idiosyncratic effects associated with each patient, either additive or multiplicative. Gaussian revival models are particularly attractive for continuous healthmeasurements because such effects are easily accommodated with block factors forpatients and temporal covariance functions such as those included in the simulationin Fig. 2.
Thus the second stage involves mainly the estimation of variance components
and range parameters in an additive Gaussian model. One slight complicationis that the revival process is not expected to be stationary, which is a relevant
consideration in the selection of covariance functions likely to be useful. Anothercomplication is that the health status may be vector-valued, Y (t) ∈ <q , so there
are also covariance component matrices to be estimated. If the covariance function
is separable, i.e.
cov(Zir (s), Zir0 (s0 )) = Σr,r0 K(s, s0 )for some q × q matrix Σ, maximum-likelihood estimation is straightforward. But
separability is a strong assumption implying that temporal correlations for all healthvariables have the same pattern, including the same decay rate, which may not be
an adequate approximation. Nevertheless, this may be a reasonable starting point.
The second stage requires all health records to be aligned at their termini. Accordingly, a record that is right censored (Ti > ci ) cannot be properly aligned. If thecomplete records are sufficiently numerous, the simplest option is to ignore censoredrecords in the second stage, on the grounds that the estimating equations based oncomplete records remain unbiased. This conclusion follows from the fact that thesecond factor is the conditional distribution given survival time. Thus, provided
that the censoring mechanism is a selection based on patient survival time, the estimating equations derived from complete records are unbiased. The inclusion ofcensored records is thus more a matter of statistical efficiency than bias, and theinformation gained from incomplete records may be disappointing in view of theadditional effort required.

5.2

Incomplete records

If we choose to include in the likelihood the record for a patient censored at c > 0, we
need the joint probability of the event T > c, the density of the subset tc = t ∩ [0, c],and the outcome Y [tc ] at y. On the assumption that censoring is uninformative,
i.e. that the distribution of the subsequent survival time for a patient censored at17

time c is the same as the conditional distribution given T > c for an uncensoredpatient, the joint density isZf (t) p(tc | t) g(y; t − tc ) dtt≥con the space of finite-length records. Assumption (3) implies that the second factor,
the density of the appointment dates in [0, c] for a patient surviving to time t > c,does not depend on the subsequent survival time t − c, in which case it may be
extracted from the integral. It is also reasonable to assume that the distributionof appointment schedules is known, for example if appointments are scheduled administratively at regular intervals, in which case the second factor may also bediscarded from the likelihood. Since the survival probability 1 − F (c) is included in
the first-stage likelihood, the additional factor needed in the analysis of the revival
model isZ
1
f (t) g(y; t − tc ) dt,
1 − F (c) t>c
in which tc may regarded as a fixed subset of [0, c]. Unfortunately, the integral
involves both the survival density f (t) = F 0 (t) and the density of the revival process,so the full likelihood no longer factors. For an approximate solution, f may bereplaced with the estimate obtained from the first-stage analysis of survival times,
and if fˆ is purely atomic, the integral is converted to a finite sum.
The situation is considerably more complicated if, as in section 4.3, the revival
processes for distinct patients are not independent,

5.3

Treatment effect: definition and estimation

We consider here only the simplest sort of revival model for the effect of treatment
on patient health, ignoring entirely its effect on survival time. Health status in therevival process is assumed to be Gaussian, independent for distinct patients, andthe treatment is assumed to have an effect only on the mean of the process, not on
its variance or covariance. Consider two patients, one in each treatment arm,ai (t) = āi (Ti − t) = 1,

aj (t) = āj (Tj − t) = 0

such that xi = xj . The revival assumption asserts that the random variable Zi (s) −Zj (s) is distributed independently of the pair Ti , Tj . By definition, the treatment
effect as defined by the revival model is the difference of meansτ10 (s) = E(Zi (s)) − E(Zj (s)) = E(Yi (Ti − s)) − E(Yj (Tj − s))at revival time s. This is not directly comparable with either of the the conventional
definitions
γ10 (t) = E(Yi (t) − E(Yj (t))

0or γ10(t) = E(Yi (t) − E(Yj (t) | Ti , Tj > t)

in which the distributions are compared at a fixed time following recruitment. The
expectation in a survival study—that healthy individuals tend to live longer than18

the frail—implies that E(Y (t) | T ) must depend on the time remaining to failure.
0
In that case, the conventional treatment definition γ10(t) depends explicitly on the
difference between the two survival times. In other words, it does not disentangle
the effect of treatment on patient health from its effect on survival time.
If the revival assumption fails, i.e. if Z is not independent of T , then, in thesimplest setting where the dependence on T is linear and additive, the difference ofmeans
E(Zi (s) | T ) − E(Zj (s) | T ) = τ10 (s) + γ(Ti − Tj ).
contains both a treatment effect and an effect due to the difference in survival times.
In other words, the failure of the revival assumption does not necessarily complicate
the interpretation of treatment effects. By contrast with standard practice in theanalysis of randomized trials with longitudinal responses, (Fitzmaurice, Laird andWare 2011, section 5.6), it is most unnatural in this setting to work with the conditional distribution given the baseline outcomes Yi (0) ≡ Zi (Ti ). That is one reason
for the recommendation in section 3.4 that the baseline response be regarded as anintegral part of the outcome sequence, not as a covariate. Exchangeability implies
distributional equality Zi (Ti ) ∼ Zj (Tj ) for individuals having the same covariate
values, but it does not imply equality of conditional distributions given T . On thepresumption that treatment assignment is independent of baseline response values,
we also have Zi (Ti ) ∼ Zj (Tj ) conditionally on treatment, whether or not ai , aj areequal. Consequently, in order to satisfy the exchangeability assumption, it is necessary to introduce a null, pre-randomization, treatment level, ai (0) = aj (0), common
to all subjects.

5.4

Revival review

The easiest way to check the revival assumption is to formulate and fit a specificalternative model in which the revival process is not independent of the survival
time. We consider here only the simplest design in which all records are complete,there are no covariates or treatment assignment, observations for distinct patientsare independent, and the revival model is a family of Gaussian process. One way
to do this is to replace (5) with
E(Zi (s) | T ) = µ(s, Ti )for some suitable family of functions µ(s, T ), leaving the covariances unchanged.
For example, if x denotes patient age at recruitment, the revival mean might bemodeled as
E(Zi (s) | T ) = µ(s) + β1 xi + β2 Tidepending additively on patient age and survival time. If β1 = β2 , the dependence
is on age at failure rather than age at recruitment. More general models involving
multiplicative interactions between s and Ti may also be considered.
Consider, for instance, the non-linear Gaussian revival model with meanµ(s) = α + βs/(γ + s),

19which is such that µ(0) = α, µ(∞) = α + β, and µ(γ) = 21 (µ(0) + µ(∞)), so that
γ > 0 is the semi-revival time. Within this family, the revival trajectory for one
patient could be different from that of another, depending on their survival times. In
other words, α, β, γ could depend on T or x + T , either of which is a violation of therevival assumption. One of the simplest models of this type is the time-accelerated
revival model in which the semi-revival time is inversely related to survival,µ(s, T ) = µ0 (sT ) = α + βsT /(γ + sT ).
As a practical matter, it would be more effective to replace γ with γ0 + γ1 /Tor exp(γ0 + γ1 /T ) to generate a test of the revival assumption. Likewise, we could
replace α with α0 +α1 Ti , asserting that the outcome sequences for long-lived patientsare elevated by a constant amount at all revival times. Similarly, if β is replacedwith β0 +β1 Ti , the the asymptote is elevated in proportion to the additional lifetime.
Any modification of this sort is a violation of the revival assumption, so thesurvival time and the revival process are no longer independent. However, thefactorization of the likelihood function remains intact, so the analysis remains relatively straightforward. For example, a likelihood ratio statistic to test the revival
assumption can be constructed by fitting two nested models to the revival process,one satisfying the revival assumption, the other involving T .

6

A worked example: cirrhosis study

6.1

Prednizone and prothrombin levels

In the period 1962–1969, 532 patients in Copenhagen hospitals with histologically verified liver cirrhosis were randomly assigned to two treatment arms, control and prednizone. Only 488 patients for whom the initial biopsy could be
reevaluated using more restrictive criteria were retained, yielding 251 and 237 patients in the prednizone and placebo groups respectively. Variables recorded atentry include sex, age, and several histological classifications of the liver biopsy.
Clinical variables were also collected, including information on alcohol consumption, nutritional status, bleeding, and degree of ascites. However, these covariates
were not included in the dataset used here, which was downloaded from the R library http://cran.r-project.org/web/packages/joineR maintained by Philipson Sousa, Diggle, Williamson, Kolamunnage-Dona and Henderson. At the end of
the study period, the mortality rate was 292/488, or approximately 60%.
The focus here is on the prothrombin index, a composite blood coagulation index related to liver function, measured initially at three-month intervals and subsequently at roughly twelve-month intervals. The individual prothrombin trajectories
are highly variable, both in forward and in reverse time, which tends to obscurepatterns and trends. In Figure 4a the mean trajectory is plotted against time fromrecruitment for four patient groups placebo/prednizone and censored/not censored,with the solid lines denoting censored individuals. Naturally, only those patients
who are still alive are included in the average for that time. Figure 4b shows thesame plots in reverse alignment. While there are certain similarities in the two plots,
20Line − Censored, Dotted Line − Failure

110(Smoothed) Mean ProThrombin Index vs. Time Until Failure/CensorshipLine − Censored, Dotted Line − Failure

120(Smoothed) Mean ProThrombin Index vs. Time

100

●
●
●
●

●

●

● ●
●

●
● ●

●

●

● ●

● ●
●
● ● ● ●
● ● ● ● ●
● ●
● ●
● ●
●
●
●
●

● ●

●

●

●

70

●
●
●

●

●
●

●

●

●

●

●
●

●
●

●

●
●

●
●

●

●
●

●
●
●

●

●

●

●

●

●

●
●

●
●

●

●
●

●

●

●

●

●

●
●

●

●

●
●

●

●

●

●
●
●

●

●

● ●

●

●

●
●

●

●

●

●

●

●

●

●

●

●

●
●

●
●

●
● ●
● ●

●

●

●

●

●

●

● ●

●
●

●

●

●
●

●

●

●

●

●

●

●

●
●

●

●
●

●

●

● ●

● ●

● ●

●
●
●

●

●

●
●

●

●

●

●

●

● ●
●
●

●

●
●
●

●
● ●

●

●

●
●
● ●

● ●
● ●
●
●

●

●

●

●

●
●

●

●
●
●

●

60

90

●
● ●

●

80

ProThrombin Index

●

● ●

●
●

●

100

●
●

● ●

●

● ●

PlaceboPrednisone

●

●

● ●
●

70
80
90
ProThrombin Index

110●

PlaceboPrednisone

50

●
●

60

●

0.0

1.0

2.0

3.0

4.0

5.0

6.0

7.0

8.0

9.0

−5

Time (in years)

−4

−3

−2

−1

0

Time Until Failure/Censorship(in years)

Figure 4: Prothrombin mean trajectories aligned by recruitment and by failure
the differences in temporal trends are rather striking. In particular, prothrombin
levels in the six months prior to censoring are fairly stable, which is in marked
contrast with levels in the six months prior to failure, as seen in the lower pair ofcurves.
Inspection of the graphs for uncensored patients in the right panel of Fig. 4, suggests beginning with the simplest revival model in which the sequences for distinct
patients are independent Gaussian with moments
E(Zi (s) | T ) = α + τāi (s) + β0 Ti + β1 s + β2 log(s + δ)cov(Zi (s), Zj (s0 ) | T ) = σ12 δij K1 (s, s0 ) + σ22 δij + σ32 δij δss0 .
The non-linear dependence on s is accommodated by the inclusion of log(s + δ) in
the mean model with a temporal offset δ, which is equal to one day in all subsequentcalculations. Inclusion of the survival time Ti is suggested by the increasing trend
along the diagonals and sub-diagonals of Table 1. Since the value at recruitment isincluded as a response for each series, treatment necessarily has three levels, null,control and prednizone. The three covariance terms are associated with independent
additive processes, the second for independent and identically distributed patientspecific constants, and the third for independent and identically distributed white
noise or measurement error. The first covariance term governs the prothrombin
sequences for individual patients, which are assumed to be continuous in time withcovariance function K1 (s, s0 ) = exp(−|s − s0 |/λ) for s, s0 > 0. The temporal range
in all subsequent calculations is set at λ̂ = 1.67 years, implying an autocorrelationof 0.55 at a lag of one year. The implied one-year autocorrelation for the observed
prothrombin sequences is considerably smaller, roughly 0.30, because of the whitenoise measurement term.
For the likelihood calculations that follow, incomplete records are ignored; only
the 1634 measurements for the 292 non-censored patients are used. The fitted
variance components, estimated by maximizing the residual likelihood, are(σ̂12 , σ̂22 , σ̂32 ) = (211.4, 209.2, 179.7),
21

all significantly positive. Using these values to determine the covariance matrix, theweighted least-squares coefficients in the mean model are shown in Table 3. The
standard error for the prednizone/control contrast is 1.77, somewhat larger than
the standard error for the prednizone/null contrast because the former is a contrast
between patients involving all three variance components, whereas the latter is acontrast within patients, which is unaffected by the second variance component.
Table 3: Regression coefficients in a revival model
Covariate
Coefficient
S.E.
Ratio
Null treatment
0.00—
Control
2.41
1.43
1.7
Prednizone
13.56
1.47
9.2
Survival (T )
1.75
0.47
3.7
Revival (s)
−2.12
0.47
−4.5
log(s + δ)
4.66
0.41
11.4
Various deviations from this initial model may now be investigated. In particular, it is possible to check whether there is an interaction between treatment
and survival time, i.e., whether the treatment effect for long-term survivors is or
is not the same as the treatment effect for short-term survivors. This comparison
involves two variance-components models having different mean-value subspaces, so
the residual likelihoods are not comparable. For likelihood comparisons, the kernel
subspace must be fixed, and the natural choice is the mean-value subspace for thenull model as described by Welham and Thompson (1997) or as implemented by
Clifford and McCullagh (2006). The likelihood ratio statistic computed in this way
is 0.83 on two degrees of freedom, showing no evidence of interaction. However,there is appreciable evidence in the data that the treatment effect (prednizone versus control) decreases as t → T , i.e. as s → 0, as is evident in Fig. 4 from theconvergence of the two lines during last six months of life. The likelihood-ratio
statistic for the treat.s interaction is 3.90 on two degrees of freedom, showing little
evidence of a linear trend, but the value for the treat. log(s) interaction is 8.68,
pointing to a non-linear trend, as is apparent from the convergence of the two lines
in Fig. 4b.
We may also check the adequacy of the assumed form for the mean model by
including an additional random deviation, continuous in reverse time, with generalized non-stationary covariance function such as K0 (s, s0 ) = −| log(s+δ)−log(s0 +δ)|.
The fitted coefficient is 2.38, and the associated likelihood ratio statistic is 1.2 on
one degree of freedom, showing no significant deviations that are continuous in reverse time. Finally, we check whether the sequences for different patients exhibita characteristic pattern or trend associated with time measured from recruitment
by including the generalized Brownian-motion covariance function −|t − t0 | in thecovariance model. The fitted variance coefficient is 2.10, and the likelihood ratio statistic is 2.38 on one degree of freedom showing no significant characteristic
patterns that are continuous in time measured from recruitment.

22

6.2

Effect of prothrombin on prognosis

Over a period of 5 years and one month following recruitment, patient u had eightappointments with prothrombin values as follows:tu (days)
Yu [tu ]0 126 226 392 770 1127 1631 185549
93 122 120 110 100
7259This is in fact the record for patient 402 who was assigned to prednizone and was
subsequently censored at 2661 days. As determined on day 1855, the survival prognosis for this patient depends on preceding sequence of measurements. Relative tothe unconditional survival density for a patient on the prednizone arm, the conditional survival density at time t > max(tu ) is modified multiplicatively by a factor
proportional to the joint conditional density of the random variable Zu [t − tu ] at
the observed point yu given Tu = t and the data observed for all other patients.
For the model described in the preceding section—in which the records for distinct patients are independent—this factor is particularly simple. The conditional
distribution of Zu [t − tu ] given T = t has a mean vector µ depending linearly ont − tu and log(t − tu + δ), and a covariance matrix Σ that is constant in t. The log
density at yu is a quadratic formh(t, yu ) = const − (yu − µ)0 Σ−1 (yu − µ)/2
depending on t only through µ. This estimated factor is shown in Fig 5a for three
versions of the record in which the final prothrombin value is 59, 69 or 79.

1.0

59

0.5

0.5

0.6

59

0.4

69

0.2

0.3

0.0

69

−0.5

79

79

5.0

5.5

6.0

6.5

7.0

5.0

5.5

6.0

6.5

7.0

Figure 5: Three versions of the record for patient 402: log modification factors forthe predictive survival density (left panel) and hazard functions (right panel).
It may be helpful to express the effect of the observed prednizone record on
the conditional survival distribution through its effect on the hazard function at23

times t > max(tu ) rather than its effect on the conditional survival density. Suppose, therefore, that the unconditional survival time for a patient on the prednizone
arm, is exponential with mean 5 years, so that the unconditional hazard function isconstant. What is the conditional hazard at time t > max(tu ) given the prothrombin sequence for patient u, with no further measurements made in the interval(max(tu ), t) other than survival? The conditional hazard functions for the subsequent two-year interval 5 < t < 7 are shown in Fig. 5b for the same three versionsof the prothrombin record. It is evident from these plots that the conditional hazard for the real patient is substantially elevated following the last measurement,
but the effect is transient and does not persist for the duration of a typical interappointment interval of one year. If the final value were 79 instead of 59, the hazard
function is almost constant, initially increasing and subsequently reverting to thelong-term value, which is slightly larger than the unconditional hazard.
The preceding analysis indicates that it may be misleading to treat the observed health sequence as a time-dependent covariate in the proportional-hazards
model. At any one failure time t measured from recruitment, some of the healthmeasurements are recent and fresh, while others are likely to be up to one year
old. Figure 5b shows that stale measurements may have negligible prognostic value
compared with fresh measurements, which contradicts the basic assumption in the
proportional hazards model. The predictive revival model automatically takes into
account the time that has passed since the last appointment, so that stale values
are discounted appropriately.

6.3

Review of assumptions

The conditional independence assumption (1) does not require appointments to bescheduled administratively, nor does it forbid patient-initiated appointments. Consider two patients i, j at time s prior to failure, having similar prior appointment
schedules and similar health values. Assumption (1) states that the conditional
appointment-initiation intensity given the observed health record and subsequent
survival time does not depend on subsequent health values. In other words, conditional independence implies that patients i, j are equally likely to initiate an appointment at time s; it is also assumed implicitly that they do so independently.
The evidence presented in section 3.2, and in Liestøl and Andersen (2002) shows
clearly that the rate of patient-initiated appointments increases in the last few
months of life. It is certainly possible that patient behaviour in this instance violates the conditional independence assumption, but the evidence presented does
not directly address the matter. All in all, assumption (1) seems unavoidable andrelatively benign.
The non-informative assumption (3) is much stronger than (1). It implies thatappointments are scheduled as if the patient will live indefinitely, which is clearly
contradicted by the evidence in this example. We now examine the consequencesof failure of (3), retaining (1).
Assumption (1) implies that the sampling is non-preferential in the sense of
Diggle, Menezes and Su (2010), which means that the second factor in (2) is thesame as if the appointment dates had been fixed by design. Consequently, the24likelihood calculations in section 5 are unaffected by the failure of (3).
If the appointment for patient u on day 1855 were self-initiated in such a waythat the last factor in (2) depends on subsequent survival, it would be technicallyincorrect to omit that factor in prognosis calculations. However, if it were known
that all appointments for patient u were on schedule, the possibility of a dependence
on subsequent survival is eliminated, and the prognosis calculations for this patientis technically correct even if the behaviour of other patients violates (3).

7

Summary

The paper examines the problem of model formulation for health sequences, whose
defining characteristic is that the state space contains an absorbing value. Each
health sequence is terminated ultimately by death, which is not equivalent to random restriction or censoring because subsequent values are known. Typically, sequence length and sequence values are not independent.
The principal suggestion is that it may be more natural in some circumstances
to align health sequences by failure time than by age or by recruitment date. The
following list describes various statistical implications of realignment.
1. The health sequence is regarded as a random process in its own right, not asa time-dependent covariate governing survival.
2. To a substantial extent, the model for survival time is decoupled from therevival model for the behaviour of the health sequence in reverse time.
3. Realignment implies that value Yi (0) at recruitment must not be treated asa covariate, but as an integral part of the response sequence. If they wereavailable, values prior to recruitment could also be used.
4. The definition of a treatment effect is not the usual one because the natural
way to compare the records for two individuals is not at a fixed time following
recruitment, but at a fixed revival time. The treatment value need not be
constant in revival time.
5. The predictive value of a partial health sequence for subsequent survival
emerges naturally from the joint survival-revival distribution. In particular,
the conditional hazard given the finite sequence of earlier values is typicallynot constant during the subsequent inter-appointment period.
6. Records cannot be aligned until the patient dies, which means that the revival
process is not observable component-wise until T is known. As a result, thelikelihood analysis for incomplete records is technically more complicated.
This aspect needs further development.
7. The omission of incomplete records from the revival likelihood does not lead tobias in estimation, but it does lead to inefficiency, which could be substantial
if the majority of records are incomplete.
8. The principal assumption, that appointment dates be uninformative for subsequent survival, does not affect likelihood calculations, but it does affect25

prognosis calculations for individual patients. For that reason, it is advisable
to label all appointments as scheduled or unscheduled.

Acknowledgements
Comments by D.R. Cox, D. Farewell, R. Gibbons, N. Keiding, S.M. Stigler andfive referees on an earlier version of this paper are gratefully acknowledged. Walter
Dempsey helped with the computation in section 6 and the preparation of Figure 4.

8

References

Andersen, P.K., Hansen, L.H. and Keiding, N. (1991) Assessing the influence ofreversible disease indicators on survival. Statistics in Medicine 10, 1061-1067.
Clifford, D. and McCullagh, P. (2006). The regress function. R Newsletter 6, 6–10.
Cox, D.R. (1972) Regression models and life tables (with discussion). J. Roy.
Statist. Soc. B 34, 187–220.
Cox, D.R. and Snell, E.J. (1981) Applied Statistics. London: Chapman and Hall.
DeGruttola, V. and Tu, X.M. (1994) Modeling progression of CD-4 lymphocyte
count and its relation to survival time. Biometrics 50, 1003–1014.
Diggle, P.J., Heagerty, P., Liang, K.-Y. and Zeger, S.L. (2002) Analysis of Longitudinal Data. Oxford Science Publications: Clarendon Press.
Diggle, P.J., Farewell, D. and Henderson, R. (2007) Analysis of longitudinal datawith drop-out: objectives, assumptions and a proposal (with discussion). Applied
Statistics 56, 499–550.
Diggle, P.J., Sousa, I. and Chetwynd, A. (2008) Joint modeling of repeated measurements and tome-to-event outcomes: The fourth Armitage lecture. Statistics
in Medicine 27, 2981–2998.
Diggle, P., Menezes, R. and Su, T-L. (2010) Geostatistical inference under preferential sampling (with discussion). Appl. Statist. 59, 191–232.
Farewell, D. and Henderson, R. (2010) Longitudinal perspectives on event history
analysis. Lifetime Data Analysis 6, 102–117.
Faucett, C.L. and Thomas, D.C. (1996) Simultaneously modeling censored survival
data and repeatedly measured covariates: a Gibbs sampling approach. Statistics
in Medicine 15, 1663–1685.
Fieuws, S., Verbeke, G., Maes, B. and Vanrenterghem (2008) Predicting renal graft
failure using multivariate longitudinal profiles. Biostatistics 9, 419–431.
Fitzmaurice, G.M., Laird, N.M. and Ware, J.H. (2011) Applied Longitudinal Data
Analysis, 2nd edition. New York: Wiley.
Fitzmaurice, G., Davidian, M., Verbeke, G. and Molenberghs, G. (2009) Longitudinal Data Analysis Chapman & Hall.
26

Guo, X. and Carlin, B. (2004) Separate and joint modeling of longitudinal and event
time data using standard computer packages. American Statistician 58, 1–10.
Henderson, R., Diggle, P. and Dobson, A. (2000) Joint modeling of longitudinalmeasurements and event time data. Biostatistics 1, 465–480.
Kurland, B.F., Johnson, L.L., Egleston, B.L. and Diehr, P.H. (2009) Longitudinal data with follow-up truncated by death: match the analysis method to theresearch aims. Statistical Science 24, 211-222.
Lagakos, S.W. (1976) A stochastic model for censored-survival data in the presenceof an auxiliary variable. Biometrics 32, 551-559.
Lee, J.C. (1988) Prediction and estimation of growth curves with special covariancestructures. J. Amer. Statist. Assoc. 83, 432–440.
Lee, J.C. (1991) Tests and model selection for the general growth curve model.
Biometrics 47, 147–159.
Laird, N. (1996) Longitudinal panel data: an overview of current methodology.
In Time Series Models in Econometrics, Finance and Other Fields, D.R. Cox,
D.V. Hinkley and O.E. Barndorff-Nielsen, eds. Chapman & Hall Monographs on
Statistics and Applied Probability 65.
Liestøl, K, and Andersen, P.K. (2002) Updating of covariates and choice of time
origin in survival analysis: problems with vaguely defined disease states. Statistics
in Medicine 21, 3701–3714.
McCullagh, P. (2008). Sampling bias and logistic models (with discussion). J. Roy.
Statist. Soc. B 70, 643–677.
Murphy, S.A. (2003) Optimal dynamic treatment regimes (with discussion). J. Roy.
Statist. Soc. B 25, 331–366.
Rizopoulos, D. (2010) JM: An R package for the joint modeling of longitudinal andtime-to-event data. Journal of Statistical Software 35, 1–33.
Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data Chapman and Hall.
Rosthøj, S., Keiding, N. and Schmiegelow, N. (2012) Estimation of dynamic treatment strategies for maintenance therapy of children with acute lymphoblastic
leukaemia: an application of history-adjusted marginal structural models. Statistics in Medicine 31, 470–488.
Schaubel, D.E. and Zhang, M. (2010) Estimating treatment effects on the marginalrecurrent event mean in the presence of a terminating event. Lifetime Data
Analysis 16, 451–477.
Sweeting, M.J. and Thompson, S.G. (2011) Joint modeling of longitudinal and timeto-event data with application to predicting abdominal aortic aneurysm growth
and rupture. Biometrical Journal 53, 750–763.

27

Tsiatis, A.A., DeGruttola, V. and Wulfsohn, M.S. (1995) Modeling the relationshipof survival to longitudinal data measured with error: applications to survival andCD4 counts in patients with AIDS. J. Amer. Statist. Assoc. 90, 27–37.
Tsiatis, A.A, and Davidian, M. (2004) Joint modeling of longitudinal and time-toevent data: an overview. Statistica Sinica 14, 809–834.
van Houwelingen, H.C. and Putter, H. (2012) Dynamic Prediction in Clinical Survival Analysis. Monographs on Statistics and Applied Probability 123; CRC
Press.
Welham, S.J. and Thompson, R. (1997) Likelihood ratio tests for fixed model terms
using residual maximum likelihood. J. Roy. Statist. Soc. B 59, 701–714.
Wulfsohn, M.S. and Tsiatis, A.A. (1997) A joint model for survival and longitudinal
data measured with error. Biometrics 53, 330–339.
Xu, J. and Zeger, S.L. (2001) Joint analysis of longitudinal data comprising repeatedmeasures and times to events. Applied Statistics 50, 375–387.
Zeger, S.L. and Liang, K.-Y. (1986) Longitudinal data analysis for discrete andcontinuous outcomes. Biometrics 42, 121–130.
Zeger, S.L., Liang, K.-Y. and Albert, P. (1988) Models for longitudinal data: ageneralized estimating equation approach. Biometrics 44, 1049–1060.

28