STATISTICS IN MEDICINE, VOL. 11,455474 (1992)

GRAPHICAL REPRESENTATION OF SURVIVAL CURVES ASSOCIATED WITH A BINARY NON-REVERSIBLE TIME DEPENDENT COVARIATE ERIC J. FEUER EPN 313, Division of Cancer Prevention and Control, National Cancer Institute, Bethesda, MD 20892. U.S.A.

BENJAMIN F. HANKEY EPN 343, Division of Cancer Prevention and Control, National Cancer Institute, Bethesda, MD 20892, U.S.A.

JEFFREY J. GAYNOR Division of Biostatistics, Memorial-Sloa Kettering Cancer Center, I275 York Avenue, New York, N . Y . 10021, U.S.A.

MARGARET N. WESLEY IMS Inc., 12501 Prosperity Drive, Suite 200. Silver Spring, MD 20910. U.S.A.

STUART G. BAKER EPN 344, Division of Cancer Prevention and Control, National Cancer Institute, Bethesda, MD 20892, U.S.A.

AND JASON S. MEYER IMS Inc.. 12501 Prosperity Drive, Suite 200, Silver Spring, MD 20910, U.S.A.

SUMMARY The use of time dependent covariates has allowed for incorporation into analysis of survival data intervening events that are binary and non-reversible (for example, heart transplant, initial response to chemotherapy). We can represent this type of intervening event as a three-state stochastic process with a starting state (S),an which usually represents death. In this paper we present intervening state (I),and an absorbing state (D), three procedures for calculating survivorship functions which attempt to display the prognostic significance of the time dependent covariate. The first method compares survival from baseline for the two possible paths through the stochastic process; the second method compares overall survival to survival with state I removed from the process; and, the third method compares survival for those already in state I at a landmark time x to those in state S at time x who will never enter state I. We develop discrete hazard estimates for the survival curves associated with the three methods. Two examples illustrate how these methods can yield different results and in which situations one might employ each of the three methods. Extensions to applications with reversible binary time dependent covariates and models with both baseline and time dependent covariates are suggested.

1. INTRODUCTION

Appropriate statistical methodology for modelling the effects of an intervening event on the hazard through a covariate whose values change after baseline, that is, a time-dependent 0277-67 15/92/040455-20$10.00 0 1992 by John Wiley & Sons, Ltd.

Received June 1990 Revised June 1991

456

E. J. FEUER ET AL.

covariate, has been in use for some time.’ The special type of time dependent covariate of interest in this paper is defined as Zj(f) =

1 for t

> wj

0 otherwise

where w jis the time subsequent to entry into a study at which individualj experiences an event or transition ofinterest. At entry into the study ( t = 0),all zj(0)= 0 and once the switch has occurred it is non-reversible. Examples of this type of time dependent covariate include intervening events that potentially decrease the hazard of death (for example, a heart transplant for individuals in a suitably defined cohort,2 the analysis of survival by tumour response3) or that increase the hazard of death (for example, the onset nephropathy among insulin dependent diabetes patients4). In the analysis of the prognostic significance of baseline covariates it is possible to obtain survivorship functions that portray the prognostic significance of the variable studied. Obtaining survivorship functions that relate to the effect of a binary non-reversible time dependent covariate is not as straightforward. Some have attempted to circumvent this problem by treating the time dependent covariate as a baseline covariate for graphical and/or analysis purposes. For example, Anderson et and others6-8 have criticized an analysis of survival by tumour response frequently used in the past in which patients are separated at baseline into two groups based on whether or not they ever achieve a response. The analysis and the associated graphical representation have an improved survival bias for the responders because patients must survive long enough to gain classification as responders. A second approach, suggested in Anderson et ~ l . is, ~to graphically compare survival for nonresponders (estimated by censoring responders at their time of response) to overall survival ignoring the intervening event. The rationale behind this approach is that when the time dependent covariate has no effect on the hazard of death these curves are equal. A final approach suggested in Anderson et al. is termed the ‘landmark method’. In this method one classifies patients at some ‘landmark time’ on the basis of whether they have responded by that time. One calculates survival conditional from the landmark time for responders and nonresponders based on their classifications at the landmark time. The purpose of this paper is to formalize and interpret these three heuristic graphical methods that represent the effect on survival of a single binary non-reversible time dependent covariate. Given the model specified in the next section of the paper, these adaptations: ~

1

.

~

3

~

1. correctly account for person-time at risk during the entire course of follow-up; 2. provide methods to calculate the hazards associated with transitions between various states; 3. show how to combine these hazards to obtain various survival curves of interest. Section 2 presents the basic data analysis problem in terms of a Cox proportional hazards model, and three survival representations associated with the corresponding heuristic methods. Section 3 presents discrete hazard estimation procedures for the three methods and demonstrates their use on the Stanford Heart Transplant data.2 Section 4 presents the graphical results of the Stanford Heart Transplant data, two hypothetical examples and data from a study of skin window reactivity to autologous breast cancer’ to demonstrate the differences in interpretation between methods 1, 2 and 3. A final section summarizes the types of questions addressed by each of the three methods and also presents some possibilities for extensions to other types of time dependent covariates.

GRAPHICAL REPRESENTATION OF SURVIVAL CURVES

457

Figure 1. The S I D stochastic process

2. THEORY 2.1. The problem

Consider a three-state stochastic process with a starting state S, an intervening state I, and an Everyone starts at time 0 in state S , and then absorbing state that usually represents death (D). either goes directly from S to D or from S to D through I. The crude hazards associated with the three possible transitions are h,,(t), h,,(t) and h,,(t), where t is the time from baseline (Figure 1). Note that we assume that the hazard associated with the transition from I to D depends only on the time from baseline and is not conditioned on the time at which the transition from S to I occurred. If we have an interest in the effect of entering the intervening state I (for example, heart transplantation) on survival, we might consider an analysis of the data that arose from this stochastic process with the use of a Cox proportional hazards model. Using this model h(tlzj(t))= ho(t)ex~ Czj(t)PI,

(1)

where h(tlzj(t))is the hazard of death for individual j at time t and Zj(t)=

i

1 if personj enters state I by time t 0 otherwise.

In the Cox proportional hazards model, a test of the null hypothesis B = 0 is equivalent to the test H,: h,,(t) = h,,(t) in the three-state stochastic model framework if we assume that h,,(t)/h,,(t) = e@for all t. To represent graphically the results of this analysis one might produce a cumulative hazard plot as in Gaynor.", l 1 Clinicians, however, generally find survival plots preferable to relate to actual patient experience. In the following sections we present three methods to graphically represent the effect of the time dependent covariate on survival. The first method represents survival through the SID stochastic process conditioned on a specific pathway (that is, the intervening event either occurs or does not occur). The second method compares overall survival through the SID stochastic process (that is, through either path) to survival through the SD stochastic process (that is, with state Z eliminated). The third method compares survival from S to D in the SD stochastic process with survival from I to D in the SZD process. Since no one is in state I at time 0, we condition the survival curves at a time x after people have entered state I to allow for the estimation of h,,(t). Simon and MakuchI2 have previously proposed this final method. We compare the three methods schematically in Figure 2.

458

E. J. FEUER ET AL.

0

Method 1

(3-03

Method 3

I=>@

0

0

I=>@ I

Figure 2. Schematic diagram of the three graphical methods

2.2. Method 1 -Path specific survival One way to represent survival for those who do and do not pass through I is to condition survival on the S-D and S-I-D pathways through the stochastic process in Figure 1 (represented as S s - D ( t ) and S s - , - D ( t ) , respectively). The probability of survival beyond time t given that one goes directly from S to D is

and the probability of surviving beyond time t given one enters I on the way to D is SS-I-D(t)

=

[:j

~ ~ s s f , s ~ ( y ) h s , ( y ) S , ,lY)h,~(x)dYdX]/8,,q (x

where: SSf,SD(t)

= exp[ -

ji

1

(hS,(x) + hSD(x))dx

(the probability of not leaving state S by time

t),

(the probability of not leaving state I by time r given entry to state I at time t,),

:1

6s-D =

~s,,sD(x)hsD(x)dx

(the probability of going directly from S to D),and 0s-i-D =

1: 1;

~si.sD(y)hsi(y)~i,(x~Y)h,D(x)dydx

(the probability of going through I on the way to D ) .

GRAPHICAL REPRESENTATION OF SURVIVAL CURVES

459

Conditioning on the pathway is a formalization of separating people at baseline into those who do and do not ever enter state 1. The advantage of the formulation presented here is that it correctly considers persons on the S-I-D pathway at risk in state S before they enter state I. As in the heuristic method, however, prospective conditioning on an event that has not yet occurred imparts an improved survival bias to those who are on the S-I-D pathway relative to those on the S-D pathway. Note that, in general, S,-,(t) # S s - , - J t ) when h,,(t) = h,,(t), and thus this method is not useful for graphical assessment of the null hypothesis in model (1). Although this method has inherent biases for comparison of survival between the two pathways, it may have use in comparison of survival on a specific path as a function of a baseline covariate (for example, males versus females on the S-I-D pathway). A slightly different approach is to compute path specific survival without conditioning on the pathway (that is, do not divide by and &-,,). This alternative approach does not yield curves scaled for direct comparison with one another (that is, at time 0 they are equal to O s - , - D and &,,, respectively, rather than both starting at 1). 2.3. Method 2 -Overall survival versus survival in the absence of the intervening event

In method 2 we compute overall survival (regardless of the path) as

Note that the lack of dashes (that is, SID instead of S-I-D) indicates reference to overall survival through the S I D stochastic process rather than through a pathway in the stochastic process. The first part of equation (2) represents the probability of transferring to state I prior to time t and surviving beyond t, and the second part represents the probability of survival beyond time t without transition prior to t. An alternative but equivalent expression for SsrD(t) is S s - , - D ( t ) ) + ( B s - D Ss- D( t ) ) , where the first term represents survival for someone who will transfer prior to death (even if it is after time t) and the second term represents survival for someone who does not. To represent survival in the SD process (that is, survival in the absence of state I) when we only observe the SID process, we must make an independent competing risks assumption between I and D from state S. This implies that the hazard of D in state S (that is, h,,(t)) remains unchanged with elimination of state I. Thus, even though we observe hsD(t)when h,,(t) > 0, we assume it takes on the same values even if the transition to state I were impossible (that is, h,,(t) = 0). We can factor the joint SDF(Ss,.sD(C)) into Ss,(t) x SsD(t) where

and,

Under the independent competing risk assumption SsD(t) represents survival with state I eliminated (that is, overall survival in the SD stochastic process).

460

E. J. FEUER ET AL.

When we compare &(t) to SsrD(t) graphically under which implies that

H, of the Cox regression model (l),

h S D ( t ) = h,,(t),

SSID(t)

=

1:

sSD(z)sSf(z)hS,(z)sS.(t

= sSD(t)

= sSD(t)

1: [1:

I z)dz + s S D ( t ) s S I ( t )

hS,(Z)SSf(Z)dz+ SSD(t)sSf(t) hSI(z)SS,(z)dz + s S I ( t )

1

= sSD(f).

Thus, this method is useful for graphically assessing the effect of the null hypothesis of model (1). In the case of an externally induced state I (for example, heart transplantation) SsrD(t) represents the survival experience of the total programme of heart transplantation, while SsD(t) represents survival with elimination of the programme. The difference between these curves is a function of not only how effective are the hearts once they have been transplanted, but how quickly hearts become available before people die without them. For an internally or biologically induced state I (for example, tumour response to chemotherapy), SsrD(t)represents survival given the pattern of responses under the chemotherapy regimen under study, while S,,(t) may represent survival if no one had a tumour response, or it may only represent a useful gauge against which one can measure the effect of the intervening event. An advantage of this formalization of the heuristic method of Anderson et aL3 is that SsrD(t) takes into account the transitions between the three states. Overall survival that ignores the underlying structure of the model (that is, use of overall survival time with disregard of the intervening event) has been shown4 to be equal to 1

where p ( x ) = the probability of being in state S at time x given one is alive at time x = SSI,SD(X)/SSI&).

In general, equation (3), which uses a weighted average of h&) and h,,(x) as the hazard function, does not equal SsrD(t). In the Appendix we present an alternative derivation of S,,(t) and Ss,,(t) that provides another interpretation for these two quantities (in which S,,, ( t ) can be interpreted as a weighted average of survival functions) and ties this method to existing methods that adjust for baseline covariates.13*l 4

2.4. Method 3-Survival representation of the covariate effect on the hazard The survival representation of model (1) by Simon and Makuch12 compares

GRAPHICAL REPRESENTATION OF SURVIVAL CURVES

46 1

to

where one chooses x arbitrarily (but ‘early’ enough so that not too many deaths have occurred prior to x and ‘late’ enough so that enough people have entered state I to obtain stable estimates of l ~ , ~ ( tSimon )). and Makuch” suggest choice of x at the median or third quartile of the ‘time to transfer’ distribution. The survival form of the time to transfer distribution is 1 I

where eS-I

=

[r

SSl,SD(y)hSI(y)dy.

Obviously, under H , SsD(tI x) = S,,(t I x). This method, however, has the disadvantage that the choice of x is arbitrary and the survival experience prior to x is not represented. Another graphical representation of the covariate effect on the hazard are cumulative hazard plots, that is

Interpretation focuses on the slope (that is, the hazard) of the cumulative hazard curve. The slope of the cumulative hazard plot is not adversely affected by having no one at risk in state I at time 0 or by the choice of a conditioning time, because evaluation of the slope commences at the point where the risk set in state I is large enough to obtain stable estimates. The survival representation is merely a transformation of the cumulative hazard plot. We term this method the ‘survival representation of the covariate effect on the hazard’, because SsD(t I x) and S,,(t 1 x) are pure functions of hsD(y) and ItID(y) for (x < y < t),respectively, and do not include any information on the transfer density into state I . In the interpretation of a survival curve, however, the focus is not merely on the slope of the curve but on its absolute level, and the choice of the conditioning point can affect this substantially. Method 3 improves upon the ‘landmark method’ in that a person who transfers to state I after the conditioning time x would have their time at risk correctly attributed to I after the transfer time, while the landmark method would attribute it all to state S. 3. DISCRETE HAZARD ESTIMATION

One possible assumption under which we may estimate the quantities in Section 2 is the discrete hazard assumption (using Kaplan-Meier” type estimates) which states that the hazard associated with an event is zero except where observed events occur. We use an early version of the Stanford Heart Transplant data as employed by Mantel and Byar’ to illustrate the methods described here. Mantel and Byar developed a Mantel-Haenszel type statistic which is equivalent to the score test of = 0 in the Cox regression model (1). The data in Table I uses a format similar to the Mantel and Byar tally sheet used in the calculation of their test statistic. To perform these calculations, we need exact transfer and loss to follow-up times not explicitly available in the ~ losses to Mantel-Byar paper. Exact transfer times were obtained from Turnbull et ~ l . , ’ and follow-up were assigned halfway between the prior and subsequent death times. One was

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

1 2 3 4 6 8 9 10 12 14 16 17 18 19 20 21 23 25 26 28 31 32 35 36 37 39 40 41 44 45 46

68 65 61 55 54 52 50 49 48 45 44 42 40 38 37 36 34 31 30 28 27 26 25 24 21 19 19 17 15 15 15

2

1 1 1

1

1

1

2 1 1

1 2 3

2

1

2

2 1

1 2 1 1 2 1 1 1 2 1 1 2 1 1 1

1

2 2 3 1

Non- transplants

o m

O~oooO 0.04oO 0.0417 00476 O.oo00 0.1053 O.oo00 O.oo00 O.oo00 O~oo00

o m

0.0250 O~oo00 O.oo00 O.oo00 0oo00 O.oo00 O.oo00 O~oooO

O.oo00

0.0227

0.0147 0.0308 0.0492 O.oo00 0.0370 00192 0.0200 ooo00 0.0208

00870 0.0500 O~oooO O.oo00 0.1176 O.oo00 O.oo00 00667

O.oooO

0.0299 00317 0.05 17 0.0182 O.oo00 0.0196 0.oooO 0.0204 0,0426 o0222 0.0233 0.0476 0.0256 0.0263 00270 0.0556 0.0294 0.0323 0.0667 0.0357 00370 0.0385

0 2 4 7 8 6 7 7 8 10 11 11 13 14 15 16 18 19 20 22 22 23 24 24 26 27 25 25 27 26 24

Table I. Illustrative tally sheet for calculating discrete hazards*

1

2

2 1

1 2 1 1 2 1 1 1 2 1 1 2 1 1 1

1

1

o m

O.oo00 O~oooO O.oo00 O.oo00 O.oo00 O.oo00 O.oo00 0.0455 O.oo00 O.oo00 O.oo00 O.oo00 O.oo00 0.0370 O~oo00 0.0800 0.0370 0.0769

o m

O.oo00 O-oo00 O.oo00 0.0909

OW00

O.oo00

0 1250

o m

O.oo00 O.oo00

2 3 1 1

undefined

2

Transplants

P

Y

h

s

rn

h,

Q\

P

50 51 58 61 66 67 68 69 71 72 77 78 83 85 100 102 149 153 188 219 252 285 340 508 675 733 852 1033 1034-t

14 13 11 10 10 10 9 9 8 7 7 7 6 5 4 4 3 2 2 2 2 1 1 1 1 1 1 1 1 1 1

1

1

1

1

~

1

OW00 OW00

1 1

O~oooo

O~oooo

o ~ o m O~oooo o ~ o ~ O~oooo O~oooo o ~ o ~ 0~5000

25 25 26 26 25 24 25 24 24 25 22 21 21 22 22 21 20 20 19 17 16 16 15 14 9 7 6 3 2 _ _ _ ~ ___

O~oooo o ~ m m O~oooo m0.m O~oooo o ~ m mo ~ m m o ~ m

m

o ~ mo ~ m o ~ m o ~ m

0900O

O*oooO

o ~ m o ~ m

0.1429 0.1667

o ~ m o m

0.1250

0.2500 0.3333

O~oooo

0.2000

O~oooo O~oooo

O~oooo

O~oooo

0.1538 0.0909

OW00

m O~oooo o ~ m o ~ m 00100 0~0000 o ~ m 0.1111 o ~ m

o

0.0714 O.oo00 OW00

1

1

2 1

1 1 1 1

1 1

1 1 1

1

2 1 1

1

1 1 1 1

1

1 1

1

1

2 1

2

2

5 1

1

1

1

1

ooooo

0.1111 0.1429 0.1667 03333

O.oooO

0.0500 0.0526 0.0588 04IOOO 00625 00667

o m o m

0.0800 0.0455 0.0476 O.oo00 O~oooO 0.0455

O.oo00

09000 09400 0.0385 0.0385 0.0400 00417 00400 09000

1034 + . Loss to follow-up times in the transplant group are 7, 39, 74, 126, 204, 252, 508, 508, 508, 508, 508, 508, 704,943,943, 1034 + , 1034 +

* Adopted from table on page 83 of Mantel and Byar.* Accession times equal transfer times. Loss to follow-up times in the non-transplant group are 23,23,

32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

W

m

P

E


1

N , ( l ) = total number in the study at time 0, and N,(i)

= N , ( i - 1)

+ T(i - 1) - F , ( i

-

1) - D,(i - 1)

for i > 1

N,(1) = 0. Using the notation from Section 2 we estimate: &'D(f(i))

= Dl(i)/Nl(i)

&,(W= W)/(jVl(~) - Dl(iN, and 61D(f(i))

= DZ(i)/NZ(i).

We define t(0)= 0 and thus &O) = h^,,(O) = h;,(O) = 0. Since deaths are assumed to occur before the transfers where there are ties, we use the notation t - ( k ) to denote the time immediately after the deaths but before the transfers, and t +(k) to denote the time immediately after both the deaths and the transfers. Thus, iSI,SD(t

-(k)) = i S I , S D ( t + ( k

iSI,SD(t

+(k)) = iSI,SD(t-(k))(l

iSI,SD(t-(o))

-

= iSI.SD(t+(o)) =

i= 1

- &D(t(k)), - iSI(r(k))?

GRAPHICAL REPRESENTATION OF SURVIVAL CURVES

465

Since M indexes the last distinct death or transition time,

if i I D ( t ( M ) l fj () ) goes to zero for any j < M (that is, if the longest follow-up time in state I represents a move to state D rather than a loss to follow-up). Os-D = 1 if S'sl,sD(r(i)) goes to zero at the longest follow-up time (that is, the longest follow-u_ptime in state S represents a transition to D or to I rather than a loss to follow-up). If i3s-D and Os-r do not sum to one then, Dinse and Larson" suggest using the modified estimates = t?-l

+

and, I

0s-D = &D/(gs-I

+ is-,).

We can interpret the 6's as the probability of going to I or D in the period of follow-up, with the @s interpreted as the probability of going to I or D eventually given that the patterns observed in the data continue past the end of follow-up. Finally, for k M ,

-=

We use the

6s here because they provide normalization constants which force

-

-

s s - D ( t ( O ) ) = ss-,(f(o)) = s S - I - D ( l ( 0 ) )

= 1.

Since the 6 s are used, these survival curves are conditioned on going on the S-D, S-I, or S-I-D pathways in the period of follow-up, and thus we define

iS-&(M)) = i s - I ( t ( M ) )= $ S - , - D ( t ( M ) ) = 0. Readers may obtain a PC program to perform these calculations written in C for IBM compatible machines by writing to the first author (EJF). 4. EXAMPLES

4.1. Stanford heart transplant study Figure 3 graphically represents survival for the data in Table I using methods 1, 2, and 3. For method 3, we use the median transfer time (25 days) as the conditioning time. The score test of

466

E.J. FEUER ET AL. Method 1

00

!

1

L

n

0

vm

in

200

Tinu (Days)

Method 2

._

1

01

Method 3 (Conditioned at Tirne=25 Days) ~ _ _ _ ~ _ _ _ _ _ _ ___

I,

02

00

1

.

0

I

n

im

n m (Days)

1

!YI

xo

Figure 3. Stanford Heart Transplant data

B = 0 in model (1) yields a p-value of 0.98. Owing to the fact that there is a waiting period for heart transplants, method 1 shows a large survival advantage for those who get a transplant. Methods 2 and 3 show survival to be approximately equal in both groups.

4.2. Hypothetical examples Consider an example where S,,(t) is exponential with the probability of death in a time interval of length 1 equal 0.1 (that is, h,,(t) = 0.105), S , , ( t ) is exponential with the probability of death in a time interval of length 1 equal 0.5 (that is, h,,(t) = 0.693), and thus h,,(t)/h,,(t) = 0.15 (proportional hazards). Furthermore, let S,,(t) be log-normal with parameters p = 3 and r~ = 0.1 and be generated independently of S,,(t). This reflects an initially increasing hazard rate that reaches a maximum value, and then decreases over time towards zero. The three survival curves appear in Figure 4.

467

GRAPHICAL REPRESENTATION OF SURVIVAL CURVES

0

10

10

30

40

Tlmo (Dayr)

Figure 4. Hypothetical Example 1: underlying survival curves

Note that the log-normal distribution with these parameters has most of the transfer time density between time 16 and 26. Figure 5 shows the survival curves associated with methods 1, 2, and 3. The median transfer time (day 20, that is, Ss-,(20) = 0 5 ) is used as the conditioning time for method 3. For this example the relative values of SsD(tlx) and SrD(tlx) do not vary with the choice of x (the conditioning time), because the assumption of exponential survival dictates that SsD(t + dlx + d ) = SsD(t I x). The probability of transferring to state I before absorption in state D,that is &-,, is 0.12. The results of the three methods differ dramatically. In method 1 conditioning on the pathway S-I-D confers a large survival advantage because there is for all practical purposes a ‘guarantee time’ of at least 16 days in state S before transferring to I. Deaths occur much more rapidly in the S-I-D pathway than on the S-D pathway after time 16, but by the time S s - , - D drops to the level of S S - D at time 29 both curves have essentially reached 0. Thus, despite the fact that being in state I confers a distinct survival disadvantage, the long waiting time in state S to get to state I provides an overwhelming survival advantage. The results of method 2 show both curves virtually identical until about time 20 after which SsrD(t) has a slight survival disadvantage relative to ssD(t), although the survival probabilities are small at this point in both curves. SsID(t)represents total survival through the SID stochastic process; thus at any time t, &(t) is a weighted average of S,&) and SrD(t)with the weights equal to the proportion in state S and I at time t , respectively (see the alternative derivation in the appendix to clarify further this interpretation of SsrD(t)).Prior to time 16 the probability of being in state 1 is essentially 0, and later, when the probability of being in state I becomes large, survival is already quite low. Thus since SsrD(t) is heavily weighted by SsD(t),SsrD(t)and SsD(t)are virtually identical except in their tails. We can conclude that survival in the SD and SID processes are similar. Despite the fact that survival in state 1 is much worse than that in state S, very few people live long enough to enter state I. 1

468

E. J. FEUER ET AL.

0

or O I

82

60

~~~

~

-

-

.

Method 3 represents the 'pure' comparison of the hazards hsD(t)and h,,(t) uncontaminated by the time to transition distribution. Thus, SsD(t(x) shows a large survival advantage over S , D ( t I ~ ) . To demonstrate how method 3 can vary as a function of the choice of the conditioning time (x) we present the following example: 0.105 0901

for t < 19 for t 2 19,

0693

for t < 19 for t 2 19,

{

hSD(t)

=

hID(t)

= 0007

i

and S,,(t) is log-normal with parameters u = 3 and o = 0-3 and is generated independently of ssD(t). As with the last example, h s D ( t ) / h [ D ( t ) = 015. The transfer distribution in the absence of D, that is, Ss,(t), has almost all of its density concentrated between times 10 and 40. The actual time to transfer distribution in the presence of D, that is, S S T r ( t ) ,has its first quartile at time 15, its

GRAPHICAL REPRESENTATION O F SURVIVAL CURVES

469

Method 3 (Conditioned at Time=15 Days)

1

h

0.8c

f

e

0.4-

0.2:

0

r , . . ...,..,,....,..,,,,...,,.,,.,..,,,......,,.,,r,J 0

0

10

20

30 T i m (Days)

! , , , , , . . . , , , , , , , , , , 1 , , , , , . . , , , , , 1 , , , , , , , , , , , , , , , , , ,

0

10

20

30

n m (Dmya)

Figure 6. Hypothetical Example 2: Method 3 Conditioned at two different times

median at time 19, and its third quartile at time 24. Figure 6 shows method 3 conditioned at time 15 and 19. Conditioning on time 15 shows a large survival advantage for not transferring, while conditioning at time 19 shows only a very modest survival advantage. Although this is an extreme example, effects of this type (but perhaps of lesser magnitude) can occur, and one should investigate the effect of the conditioning time on the resultant survival curves when using method 3. 4.3. Skin window study

Black et aL9 studied the relationship of cell-mediated immunity to autologous breast cancer tissue and subsequent clinical behaviour. The cell mediated immunity is measured by the cellular response to coverslip-mounted sections of autologous breast cancer tissue applied to a microabrasion of the skin (that is, a ‘skin window’ (SW)). The tests were conducted at irregular intervals after breast cancer surgery and, at the time of the first test, patients had to be free of any evidence of recurrent disease or second primary breast cancer. They showed that a positive SW reading (SW + ) is a significantly favourable prognostic factor for systemic recurrence-free survival. The data were first analysed by treating the time dependent covariate as a baseline covariate (that is, a patient was considered SW + from baseline onward if any test in the first two years was

470

E. J. FEUER ET AL.

Method 2

I

i

(0

01

08

0.

02

o*

-

Method 3 (Conditioned at Tim=5 Months)

I

-

-

.

lm

0

Tim. (hnthm)

Figure 7. Skin window data

positive and SW - if all tests in the first two years were negative) so that one could ‘graphically portray the biological significance of SW reactivity’ (p. 76). This representation has a biased survival advantage for SW patients because they must remain systemically recurrence free long enough to become SW + . Later in the paper, there was an analysis of the data using the time dependent covariate

+

Zj(t) =

i

1 if any SW test for womanj in the interval (0, t ] is positive 0 if all SW tests for womanj in the interval (0, t ] are negative,

although with no associated graphical representation. Figure 7 represents the three methods applied to the SW data where state S corresponds to a woman who is SW negative on the first post-operative test (within 6 months of surgery), a woman enters state I at her first SW positive reading, and D represents the occurrence of systemic metastases or death from disease without prior detection of systematic metastases. If the first test (within 6 months of surgery) is positive,

GRAPHICAL REPRESENTATION OF SURVIVAL CURVES

47 I

the woman transfers to the SW positive state immediately after surgery. With method 3 we used the median transfer time (5 months) for the conditioning time. This analysis is similar to the time dependent covariate approach taken by Black et While Black et al. assumed that a woman is not in either risk set (SW + or SW -) until she has her first test, we wanted to classify a woman immediately after surgery. To do this we: 1. excluded women who did not have a skin window test within six months of surgery for breast cancer (this reduced the sample size from 267 to 21 1); 2. assumed that the skin window status at the first test after surgery applied from the time of surgery.

Despite these differences, and the inclusion of the baseline covariates nuclear grade and node status in their model, we both obtained a hazard ratio of 0.27 associated with a survival advantage of skin window positivity. &I = 0.53 which indicates that we would expect 53 per cent to transfer to the SW + state within 136 months (the longest follow-up time in the SW - state). Figure 7 shows the graphical representation of the three methods. Method 1 shows the largest survival advantage for skin window positivity, but is again influenced by the artifactual influence of a ‘guarantee’ time until the change to a positive status. Method 2 shows that overall survival would decrease slightly if no one ever became SW + ,but the size of the difference is smaller than that observed in method 1. This is influenced by the fact that 47 per cent do not become positive within the period of follow-up, the shape of the transfer hazard from SW - to SW +, and the relative size of the h,, and h,, hazards. Method 3 isolates this last factor and confirms the fairly large survival advantage which skin window positivity confers. 5. DISCUSSION The three methods presented in this paper correctly account for person-time at risk as a person transfers between states, and the actual state the person is in when failure occurs. These methods address three different questions one can pose, and the choice of method to be used depends upon the substantive question of interest. Method 1 asks: ‘What is survival conditional on a specific path?’ Since we never know the path prospectively, this method does not yield a useful graphical representation of the effect of the time dependent covariate. We may, however, find it informative to compare path specific survival as a function of a baseline covariate (for example, males versus females on the S-I-D pathway). Method 2 asks: ‘What is the survival impact of elimination of state I from the SID process?’ It is affected by both the survival in state I once the transfer is made, as well as the probability of transfer and the shape of the transfer density. Method 3 asks: ‘What is the direct survival comparison associated with the hazards h,,(t) and h,,(t)?’ It is unaffected by the probability of transfer and the transfer density. In the heart transplant example method 2 examines the efficacy of the entire transplant programme (versus the alternative of eliminating the programme) while method 3 examines effectiveness of hearts once transplanted. Both methods 2 and 3 require satisfaction of the assumption of independent competing risks for &(t) to represent survival from S to D in the absence of I. If this assumption is not met the graphical representation remains valid, that is, in both methods under the null hypothesis of fi = 0 in model (1) the two curves will be equal, and method 3 remains a transformation of the cumulative hazard plots. Under non-independent competing risks S,,(t) still represents a survival curve that constitutes a useful standard from which to measure the effect of an intervening event. There are several possible extensions of the basic model we have developed for a binary non-reversible time dependent covariate. An immediate extension is to allow a person to start in state I at time 0 by simply letting N,(l) in Table I equal the appropriate number at risk at time

472

E. J. FEUER ET AL.

t( 1). Another possible extension is to a reversible binary time dependent covariate. The only

necessary addition to Table I is to include a column for transfers out of the transplant group and to allow for accessions into the non-transplant group. The complication comes in defining in method 1 which of the multitude of possible paths prior to death one wishes to compare. One would calculate overall survival through the entire process in method 2 as the time to absorption distribution in a semi-Markov process.’ ’-I9 Other extensions which are possible within the Mantel-Byar framework include: (1) time dependent covariates that have more than two levels, and (2) stratification according to discrete baseline covariates. Extension (1) is achieved by simply keeping track of the various transitions and person-time at risk in the various states (although the number of possible pathways through the process quickly becomes large), while extension (2) is achieved by performing a separate analysis in each strata defined by the cross-tabulation of the baseline covariates. Both are limited by the availability of data to produce stable estimates of the transition hazards. Extension of this methodology to include continuous baseline covariates requires explicit use of estimates from either a parametric model or the semi-parametric Cox model. For example, if there is a binary time dependent covariate ( z j ( t ) and ) a continuous baseline covariate ( Wj), we first fit the following two Cox regression models: hl(flzj(t)i Wj) = ho(t)exPCzj(t)B1 + WjDzl

where the event of interest is death, and

h,(t I Wj) = go(t)exPC w j ~ l where the event of interest is a transition to state I , and g O ( t )is the underlying hazard. For a fixed value of W = k we can estimate the three hazards necessary to calculate methods 1,2, and 3, that is, k D ( t I W = k ) = ho(t)expCkB21, hID(t

I

hSAt I

= k, = hO(t)exp[Bl

+kp219

w = k ) = so(t)expCkyl.

Kalbfleisch and PrenticeZo(p. 84-86 and 134) present the methodology for deriving the baseline hazard function, h,(t), in the presence of time dependent covariates. Finally, for a continuous time dependent covariate there are an infinite number of possible pathways through the process, and we suggest discretizing the covariate prior to analysis to limit the number of possible pathways. APPENDIX: AN ALTERNATIVE DERIVATION OF METHOD 2 Consider the random variables T as the total survival time, Ts, as the ‘time due to go from S to I’, and T S D as the ‘time due to go from S to D’. Both Ts, and T S D are hypothetical or potential times in that one will observe only the smaller of the two. A person follows the path S-I-D if Ts, < T S D and the path S-D if T,, > T S D . Then,

and, Ss-j-D(t)

= P ( T > t I Ts,
t,

TSD).

GRAPHICAL REPRESENTATION OF SURVIVAL CURVES

473

and SS,(t) = P ( T S I

> t).

We can also interpret ssD(t) as P ( T > t I Ts, = co )which is the SDF for death with the potential to go to I eliminated. Next, given that Ts, and TsD are independent, we can define SsrD(t I 20) = P ( T > t I Ts, = zo) =

{

[SSD(zO)SID(t 1 zO)l CsSD(t)l

if > zO if t < zo'

Note that one can die prior to z,; the fact that TsI = zo only says that one will transfer to I if one lives that long. If we weight SsrD(t1 z) by the potential transfer time density fsr(z) = - dCSs,(z)l/dz

= Ss,(z)hs,(z)

then

This latter quantity is identical to (2) in Section 2.3. This method of obtaining survival curves by taking a weighted average is similar to the 'direct adjustment' approach employed by MakuchI3 and Chang et ~ 1 . to ' ~ obtain survival curves for different treatment groups weighted by a common distribution of other possibly confounding baseline covariates. In their situation they have two (or possibly more) covariates with: 0 for a control group, z l = { 1 for a treatment group,

and z2 = i representing the levels of a confounding covariate. We represent survival for the two treatment groups by S,(t) =

1i P ( T > t I z1 = 0,z2= i ) P(z2 = i )

S 2 ( t )=

1P ( T > t 1 z1 = 1, z 2 = i ) P ( z 2 = i).

and, i

One advantage of this method of representation of SsrD(t) is that we can apply weighting functions other than fs,(z) to SsID(tI z). For example, one might want to know the effect of a concerted media effort to make more hearts available which would change the potential transplant time distribution.

474

E. J. FEUER ET AL. ACKNOWLEDGEMENTS

The authors wish to thank Dr. Maurice Black and Dr. Reinhard Zachrau, Institute of Breast Diseases, New York Medical College, Valhalla, N.Y. for the use of their data in the example of Section 4.3. Dr. Gaynor’s effort in this research was funded in part by Grant No. CA-47179 from the National Institutes of Health. REFERENCES 1. Cox, D. R. ‘Regression models and life tables (with discussion)’, Journal of the Royal Statistical Society, Series B, 34, 187-220 (1972). 2. Mantel, N. and Byar, D. P. ‘Evaluation of response - time data involving transient states: An illustration using heart transplant data’, Journal of the American Statistical Association, 69, 81-86 (1974). 3. Anderson, J. R., Cain, K. C. and Gelber, R. D. ‘Analysis of survival by tumor response’, Journal of Clinical Oncology, 1, 710-719 (1983). 4. Andersen, P. K. ‘Multistate models in survival analysis: a study of nephropathy and mortality in diabetes’, Statistics in Medicine, 7, 661-670 (1988). 5. Anderson, J. R., Cain, K. C., Gelber, R. D. and Gelman, R. S. ‘Analysis and interpretation of the comparison of survival by treatment outcome variables in cancer clinical trials’, Cancer Treatment Reports, 69, 1139-1144 (1985). 6. Weiss, G. B., Bunche, H. 111 and Hokanson, J. A. ‘Comparing survival of responders and nonresponders after treatment: a potential source of confusion in interpreting cancer clinical trials’, Controlled Clinical Trials, 4, 43-52 (1983). 7. Mantel, N. ‘Responder versus nonresponder comparisons: danuorubicin plus prednisone in treatment of acute nonlymphocytic leukemia’, Cancer Treatment Reports, 67, 3 15-316 (1983). 8. Oye, R. K. and Shapiro, M. F. ‘Reporting results from chemotherapy trials: Does response make a difference in patient survival?, Journal of the American Medical Association, 252, 2722-2725 (1984). 9. Black, M. M., Zachrau, R. E., Hankey, B. F. and Wesley, M. ‘Skin window reactivity to autologous breast cancer, an index of prognostically significant cell-mediated immunity’, Cancer, 62, 72-83 (1988). 10. Gaynor, J. J. ‘The use of time dependent covariates in modelling data from an occupation cohort study’, Journal of the Royal Statistical Society, Series C, 36, 340-351 (1987). 1 1 . Gaynor, J. J. ‘Comment on “A reanalysis of the Stanford heart transplant data” by Aitkin, Laird, and Francis’, Journal of the American Statistical Association, 78, 288-290 (1983). 12. Simon, R. and Makuch, R. W. ‘A non-parametric graphical representation of the relationship between survival and the occurrence of an event: Application to responder versus non-responder bias’, Statistics in Medicine, 3, 35-44 (1984). 13. Makuch, R. W. ‘Adjusted survival curve estimation using covariates’, Journal of Chronic Diseases, 35, 437-443 (1982). 14. Chang, I. M., Gelman. R. and Pagano, M. ‘Corrected group prognostic curves and summary statistics’, Journal of Chronic Diseases, 35, 669-674 (1982). 15. Kaplan, E. L. and Meier, P. ‘Nonparametric estimation from incomplete observations’, Journal of the American Statistical Association, 53, 457-48 1 (1958). 16. Turnbull, B. W., Brown, B. W. and Hu, M. ‘Survivorship analysis of the heart transplant data’, Journal of the American Statistical Association, 69, 74-80 (1974). 17. Dinse, G. E. and Larson, M. G. ‘A note on semi-Markov models for partially censored data’, Biometrika, 73, 379-386 (1986). 18. Lagakos, S. W., Sommer, C. J. and Zelen, M. ‘Semi-Markov models for partially censored data’, Biometrika, 65, 311-317 (1978). 19. Howard, R. A. Dynamic Probabilistic Systems Vol. II Semi-Markou and Decision Processes, Wiley and Sons, New York, 1971. 20. Kalbfleisch, J. D. and Prentice, R. L. The Statistical Analysis of Failure Time Data, Wiley, New York, 1980.