Semiparametric methods for multistate survival models in randomised trials.

Research Article Received 30 May 2012,

Accepted 18 November 2013

Published online 13 December 2013 in Wiley Online Library

(wileyonlinelibrary.com) DOI: 10.1002/sim.6060

Semiparametric methods for multistate survival models in randomised trials‡ Harold M. Hudson,a,b,c * † Serigne N. Lô,b R. John Simes,a Andrew M. Tonkind and Stephane Heritierb Transform methods have proved effective for networks describing a progression of events. In semi-Markov networks, we calculated the transform of time to a terminating event from corresponding transforms of intermediate steps. Saddlepoint inversion then provided survival and hazard functions, which integrated, and fully utilised, the network data. However, the presence of censored data introduces significant difficulties for these methods. Many participants in controlled trials commonly remain event-free at study completion, a consequence of the limited period of follow-up specified in the trial design. Transforms are not estimable using nonparametric methods in states with survival truncated by end-of-study censoring. We propose the use of parametric models specifying residual survival to next event. As a simple approach to extrapolation with competing alternative states, we imposed a proportional incidence (constant relative hazard) assumption beyond the range of study data. No proportional hazards assumptions are necessary for inferences concerning time to endpoint; indeed, estimation of survival and hazard functions can proceed in a single study arm. We demonstrate feasibility and efficiency of transform inversion in a large randomised controlled trial of cholesterol-lowering therapy, the Long-Term Intervention with Pravastatin in Ischaemic Disease study. Transform inversion integrates information available in components of multistate models: estimates of transition probabilities and empirical survival distributions. As a by-product, it provides some ability to forecast survival and hazard functions forward, beyond the time horizon of available follow-up. Functionals of survival and hazard functions provide inference, which proves sharper than that of log-rank and related methods for survival comparisons ignoring intermediate events. Copyright © 2013 John Wiley & Sons, Ltd. Keywords:

cumulative incidence; CIF extrapolation; Kaplan–Meier integral; multistate model; semi-Markov network; saddlepoint approximation

1. Introduction Multistate models (MSMs) involving networks (or flowgraphs) of intermediate (transient) states before a terminal event (defining one or more absorbing states) have considerable appeal in the analysis of randomised controlled trial (RCT) data. This is because successive states are clinically meaningful in the study of disease progression. In this paper, we provide a semiparametric alternative to our earlier parametric approach [1], which employed saddlepoint inversion to estimate survival and hazard functions of time to trial endpoint in semi-Markov networks. The parametric approach requires specifying parametric distributions of gap times to successive transitions in a multistate network. In this setting, the time from study entry to the next succeeding event (whether intermediate or final outcome) is often termed the time to first event in clinical medicine literature, whereas it is known as the initial state holding time in engineering and other contexts. For later transitions, including those with competing event outcomes, we shall refer to time to a NHMRC

Clinical Trials Centre, University of Sydney, New South Wales 2006, Australia George Institute, University of Sydney, New South Wales 2050, Australia c Department of Statistics, Macquarie University, North Ryde, New South Wales 2109, Australia d Department of Epidemiology and Preventive Medicine, Monash University, Melbourne, Victoria 3004, Australia *Correspondence to: Harold M. Hudson, Department Statistics, Macquarie University, North Ryde, New South Wales 2109, Australia. † E-mail: [email protected] ‡ SIM Home Page at www.interscience.wiley.com/ jpages/ 0277-6715 b The

Statist. Med. 2014, 33 1621–1645

1621

Copyright © 2013 John Wiley & Sons, Ltd.

H. M. HUDSON ET AL.

1622

next event. For a semi-Markov model, the time scale applying to the hazard of transition from a transient state begins from entry to that state, that is, employing a clock-reset approach [2]. Because of the difficulty of correctly specifying component distributions, there is appeal in nonparametric estimation of the holding times in states of the network. Butler and Bronson [3] have recently adapted saddlepoint methods developed in a series of earlier works [4–6] for nonparametric estimation for survival and hazard functions of time to endpoint. In combination with efficient bootstrap approaches employing saddlepoint methods, we developed confidence envelopes for estimates. However, in trials with censoring of survival times extending beyond an administrative cut-off, the nonparametric methods of Butler and Bronson do not readily allow for nonparametric survival or hazard estimation. The purpose of this paper is to describe this context, and to propose and test methods that bridge the gap between parametric and nonparametric methods suitable for such clinical trial contexts. For nonparametric estimation, incomplete survival data introduce difficulties, because censored cases cannot be assigned a subsequent pathway to the endpoint, and progression through states after the trial’s period of follow-up is unknown. For example, with type 1 censoring, a substantial technical problem arises for transform methods: the complete data empirical estimators of a generating function cannot be consistently estimated beyond the period of the trial. This is because some of the study population may remain alive at closure of the trial, so that empirical estimators of transforms of time from entry to terminal state require unavailable survival experience past a close-out time, , specified in the study design. In clinical trials, this loss to follow-up often occurs as end-of-study censoring (also known as administrative censoring, with special case type I censoring) that limits the period of observation (to maximum time on study, the valid domain of the empirical survival function). For trials with an effective maximum period of follow-up, we introduce and recommend model-based estimation of subsequent survival and cumulative incidence functions (CIFs). The first objective of this work is to extrapolate empirical estimates available in observed periods of follow-up in each component of the network so that transforms are well defined. Although the basis of a general approach is provided here, in order to provide closed-form transforms of the extrapolated survival distributions, we conveniently chose the hazard rates for residual times beyond the period of observed follow-up as constant or linear in time. This approach differs from a common method, which employs a rescaling technique [7, 8]. Because our extrapolation approach makes assumptions about survival beyond study completion, which are untestable using trial data, it is important to gauge the impact of end-of-study censoring on transform methods. Sensitivity analysis of the impact of any such model-based extrapolation on survival estimation and inference is appropriate. It is possible that survival estimates in the observed period of follow-up will be little affected by the specific extrapolation employed. If so, semiparametric transform methods will provide robust estimation and inference in network MSMs without sacrificing the fidelity to data. Otherwise, it is important to provide methods that allow the investigation of survival beyond time under various scenarios. Evaluation of the sensitivity of inference to model specification is the second objective of our paper. Cumulative incidence functions and their empirical estimators are well known in contexts with competing risks; see equations (12) and (17) of [2]. CIFs are necessary to define transforms of state occupancy times for states from which exits can occur to one of several alternative states at the time of transition. These states can be considered to be competing causes of exit. Methods of extrapolating CIFs have only recently been discussed in the literature. Vertical modelling [9] is an approach, which combines estimation of survival to next event with the relative hazard of events of each type at the event time; see particularly equation (13). Extrapolation of a CIF beyond time is then provided by specifying a parametric extrapolation of cumulative hazard of any transition together with a relative hazard of each event type at this transition. We apply constant relative hazard extrapolations in transform methods and saddlepoint inversion. The flexibility of our approach makes it suitable for sensitivity analysis of MSM findings (refer to Section 4.4). As with fully parametric models, saddlepoint inversions provide both the survival function and hazard function of first passage time to a terminal state. These estimators fully utilise the data in the network, being derived from a matrix of transmittances, which summarise information available in network pathways leading to the trial endpoint. Functionals of survival and hazard functions may be estimated within trial treatment groups, together with their resampling distribution. Such functional estimates can be used to calculate net treatment benefits within the network, and assess evidence of non-proportional hazards (non-PHs). In Section 2, we consider the approach of [3] to survival and hazard estimation for censored data networks. We described its application in a specific trial context, employing a three-state illness-death Copyright © 2013 John Wiley & Sons, Ltd.

Statist. Med. 2014, 33 1621–1645

H. M. HUDSON ET AL.

Figure 1. Network diagram of Long-Term Intervention with Pravastatin in Ischaemic Disease study events showing randomisation-stroke-death sub-model.

model. Semi-Markov processes (SMPs) differ from continuous-time Markov chains by assuming that waiting time distributions are arbitrary, not memoryless, and depend on the destination state. The independence assumptions between waiting times allow us to express transforms of overall survival time in terms of transforms of the network components. For example, in the simple illness-death model, the moment-generating function (MGF) of overall survival time is simply expressed in terms of the MGFs of time durations of steps from states 1 to 2, 2 to 3 and 1 to 3. This simplification applies to any finite-state semi-Markov network and the algebraic relationship expressed in terms of a transmittance matrix. We defined transmittance from one state to an immediate successor as the product of the probability that the path is followed and the MGF of time before exit, given the successor state. Similarly, we defined transmittance for first passage time through a network to a designated outcome as the product of the conditional probability of the outcome, given initial state, and the MGF of first passage time to the designated outcome, given that it is reached. We recovered accurate survival and hazard function estimates, together with resampled estimates providing confidence envelopes, by saddlepoint inversion of the transmittance from entry to endpoint, determined from the transmittance matrix and an algebraic relationship (Butler and Bronson, Theorem 1). In Section 3, we consider nonparametric analysis with censored observations in more detail. With type 1 censoring, novel extrapolations of hazard of first event and cause-specific event incidence introduced in this section are necessary. To estimate survival and hazard functions of time to endpoint in semi-Markov multistate networks, empirical transform estimates can replace transforms of parametric distributions. We provide an approach avoiding the assumption of tail-estimable censoring, a condition required for consistency of saddlepoint estimators, which is often not met in RCTs. Figure 1 displays a network relevant to analysis of the Long-Term Intervention with Pravastatin in Ischaemic Disease (LIPID) study [10]. This study is particularly appropriate to examine because of its very large nature, having randomised in excess of 9000 subjects who were followed for a mean of 6.0 years, with ascertainment of vital status at completion in all but one patient and adjudication of all major events, including those tested here. We gave our analysis of the LIPID study utilising information on stroke events occurring during follow-up on a pathway (similar to that shown with solid line arrows in the figure) to death of any cause in Section 4. This analysis illustrates particular benefits of MSMs and transform methods: (i) estimation of distributions of waiting time to terminal events independently, by treatment; and (ii) valid and powerful treatment comparisons without dependence on PH assumptions. In comparison, a standard Cox model of time to study endpoint ignores intermediate events in primary endpoint analysis.

2. Semi-Markov multistate networks


Statist. Med. 2014, 33 1621–1645

1623

Multistate networks are defined in clinical trials by the occurrence of intermediate and terminal events followed in the study. Occurrence of an event initiates a new health state for the individual trial participant. With this transition is associated the corresponding time duration spent in the preceding state. Such networks may contain three components of special interest: sequences of states constituting a pathway to outcome, competing events and loops. In what follows, we describe an irreversible network in

H. M. HUDSON ET AL.

which loops are absent. In such cases, return to a previous state cannot occur. Although this assumption simplifies discussion and interpretation, estimates are available using the approach to be described for networks including loops [5, 11]. A direct pathway to an outcome specifies a sequence of necessary intermediate events, whereas competing events define alternative successor states and pathways. With semi-Markov assumptions, we characterised the transmittance of time from entry to terminal event by the probabilities of passage to alternative successor states and the MGFs of distributions of time spent in the initial and intermediate states before progressing. 2.1. Semi-Markov model of multi-stage Long-Term Intervention with Pravastatin in Ischaemic Disease outcomes Consider the LIPID flowgraph (Figure 1). The figure shows alternative pathways to the endpoint of a major scientific study of a cholesterollowering treatment. The LIPID study investigated effects of randomisation to pravastatin compared with placebo on pre-specified outcomes (Table I). If the outcome being analysed is death of any cause, and intermediate acute myocardial infarction events are not considered, the LIPID network simplifies to the randomisation-stroke-death sub-model, which is a simple illness-death model. In this model, a variety of individual event histories can occur, with corresponding exits from states, as illustrated in Figure 2. For clarity in exposition, our discussion and analysis will use this sub-model. We chose our focus on all-cause death, rather than death attributed to coronary heart disease (CHD) cause, to avoid the small additional complication of considering competing risks after the intermediate event stroke. Although the findings in the LIPID study were similar for these two outcomes, all-cause death is an unequivocal endpoint. Thus, the form of event histories considered are those illustrated in Figure 2. Here, events of interest are stroke and death, with end-of-study censoring occurring soon after—although not necessarily at—a planned time . Instead, the distribution of administrative censoring time has finite support with upper limit H , whereas the longest observed period of follow-up is O . Our analysis will use times to next event (exit times from states in the pathways to endpoint). Maximum observed follow-up for exits from states 1 and 2 is denoted by O1 and O2 . Our method sets cutpoints, here 1 and 2 , for the purpose of extrapolation. In sub-figure (c), the time scale has been reset after the first event (stroke) experienced by patients 2, 3 and 4. The network diagram illustrates the reality observed in many RCTs: multiple events, sometimes recurrent, with presence of competing risks. Analysis of randomised trial data often employs the Cox model for survival time to a specified outcome. The alternative MSM employs a network of pathways, such as those in Figure 1, to the study outcome. The choice is between a single summary, normally the hazard ratio (HR), and reporting of network parameters. The ability to not only ascertain but also interpret the treatment effect, together with the opportunity to assess a randomised comparison of outcomes, is important to clinical acceptance of the findings from statistical analysis. But parameters and coefficients of MSM analyses apply to single transitions within the network. We believe that integration of network elements in MSMs is therefore a necessary requirement for acceptance of such analyses in clinical research employing RCTs. The PH assumption for time to all-cause death is not required by analysis using transform methods. Survival and hazard estimation can be conducted separately in each treatment arm, or analysis can be conducted in a single arm (or cohort). Alternatively, any standard approach to MSM parameter estimation

Table I. Long-Term Intervention with Pravastatin in Ischaemic Disease pre-defined outcomes. Event

1624

Death from CHD cause Death from any cause Death from cardiovascular cause Coronary event (CHD death and non-fatal AMI) AMI, stroke, unstable angina, coronary revascularisation

Comment Primary outcome Secondary outcome Secondary outcome Composite secondary outcome Intermediate events, which are secondary outcomes

CHD, coronary heart disease; AMI, acute myocardial infarction.


Statist. Med. 2014, 33 1621–1645

H. M. HUDSON ET AL.

Figure 2. Sample of observed data patterns for five patients; sub-figure (a) represents the pathway through the whole network (1-2-3): patient 1 was event-free for 6.25 years after entry (randomisation) at which time the trial closed (end of follow-up, ‘eof’); patients 2 and 3 experienced stroke and remained alive during follow-up; patient 4 also experienced stroke but died 6 months later; patient 5 died 2 years after entry. Sub-figures (b) and (c) show the exit time from state 1 (for all five patients) and from state 2 (for three patients who experienced stroke as first event), respectively. Times 1 and 2 denote chosen cutpoints for the time to first event and the time to death after entering state 2, whereas O1 and O2 denote the largest observed time duration before exit from states 1 and 2, respectively.

(e.g. estimation employing a covariate indicating the randomised treatment arm) can provide empirical survival estimates for the method. However, in trials with a network of events defining an MSM, it seems unlikely, a priori, that PHs of the survival to terminating event will hold. This is because a common HR is unlikely to apply in all stages. Indeed, in the LIPID trial, the ratio of cumulative hazards of death in pravastatin and placebo study arms—shown in Figure 3—appears time-dependent, consistent with accumulating benefit from continued pravastatin use. An MSM in the LIPID study analysis provides the ability to remove or relax PH assumptions. In order to better interpret the output of an MSM, questions involving transition probabilities in the network are formulated. In the LIPID study, MSMs describe each transition and corresponding hazard rates of successor events. Model fitting within the randomisation-stroke-death sub-network of Figure 1 can clarify: whether hazard of death is proportional in the treatment arms; whether any evidence of accumulating benefit of statin use is statistically significant; whether the probability of following the stroke pathway differs by treatment arm; whether improved overall survival in those randomised to a statin is attributable to reductions in (first event) stroke alone; and whether the relative incidence of stroke to death remains constant (in time or by treatment). The inclusion of time to event data for intermediate events together with the use of corresponding inferential methods enables this clarification. 2.2. Transform methods


Statist. Med. 2014, 33 1621–1645

1625

In the analysis of RCTs, explicit estimation of the survival function of time from randomisation to the study endpoint remains of primary importance when conducting an MSM fit. Saddlepoint inversion will

0.20

H. M. HUDSON ET AL.

0.05 0.01

0.02

Cum.Hazard

0.10

controls pravastatin

1

2

3

4

5

6

time on study (yrs)

Figure 3. Survival time from randomisation, by allocated treatment (pravastatin or placebo control). Proportional hazards would result in parallel curves of cumulative hazard (log scale).

provide estimates of both survival and hazard functions once the MGF of the distribution of time from randomisation to endpoint is specified. Transform methods provide this specification, in terms of MGFs of network survival time components. To introduce transform methods, restrict attention to the randomisation-stroke-death sub-network of Figure 1 with patients following two pathways following randomisation (entry to state 1) to eventual death (terminal state 3). A patient’s life history following randomisation is either direct to death or with intermediate stroke (transition to stroke, state2). To estimate survival through the system, observe that SY .t / D pS123 .t / C .1 p/ S13 .t / where S123 is the waiting time distribution from randomisation to entry into state 3, conditional on the pathway being 1!2!3 and p is the probability of this pathway being followed. Note that this example has the two structural components of irreversible networks, a node (here, node 1) with waiting time to first of two competing risks, and a pathway (1!2!3) in which a sum of survival time components is needed. With semi-Markov assumptions applying to the network, the MGF of any sum of state holding times along a pathway is simply related to the MGFs of each of its components. In the randomisation-strokedeath model, with passage from 1 ! 2 ! 3 with probability p with independent transition times with survival functions S12 and S23 , and passage from 1 to 3 with probability 1 p, the MGF of survival time is MY .s/ D pM12 .s/ M23 .s/ C .1 p/M13 .s/, where Mij represents the MGF of holding time in state i before exit to state j . The terms pM12 .s/; M23 .s/ and .1 p/M13 .s/ are termed transmittances of corresponding network branches. The transmittance from entry to endpoint is MY .s/ because the probability of death is 1. We shall use the notation Tij .s/ for the transmittance for the first passage time from state i to state j . Similar expressions are readily obtainable for any finite-state semi-Markov network by a calculation using cofactors of a transmittance matrix [6, Chapter 13]. 2.3. Saddlepoint inversion

1626

Saddlepoint inversion [3, 6, 12] uses the cumulant-generating function (CGF), log MY .s/, of time to outcome to provide accurate approximations of survival and hazard functions. For simple illness-death models, exact nonparametric maximum likelihood (NPML) estimates of the survival function exist under conditions given in [13]. Exact expressions for the survival function in SMPs are also available using the semi-Markov renewal equation. However, saddlepoint approximations are generally available and are known to be highly accurate in suitable contexts (see discussion in the succeeding text under Section 5). The approach is computationally efficient, with important benefits for complex networks and bootstrap resampling schemes. In an SMP network, saddlepoint estimation is a very convenient method of approximating not only the distribution of passage time to the specified endpoint but also its hazard function. After fitting an MSM for the network, a saddlepoint expansion refocusses attention on survival and hazard to trial endpoint, thereby integrating findings for component elements of the network. Copyright © 2013 John Wiley & Sons, Ltd.

Statist. Med. 2014, 33 1621–1645

H. M. HUDSON ET AL.

2.4. Estimation Likelihoods and estimation must account for censoring, assumed to be independent of the SMP. Maximum likelihood estimation is used with parameterised network component distributions, whereas semiparametric methods (using empirical survival and competing risks estimators) are a focus of our contribution in Section 3. The likelihood of component holding times in states of SMPs is readily formed [7]. This likelihood factors conveniently to permit estimation in each component transition of the network model. Exit times along a path to a specified destination state are sufficient for the relevant model transition parameters. The presence of a competing risk, as in first-event analysis of an illness-death model, provides a mixture distribution of waiting time Y through the system. [8] demonstrated that NPML estimation of competing risk sub-distributions, S12 and S13 , are provided by the empirical CIFs for each competing cause (in the LIPID model: death, stroke). A stratified analysis, repeating analysis separately in each trial arm, is feasible following estimation of MSM component distributions. 2.5. Effect of type 1 censoring


Statist. Med. 2014, 33 1621–1645

1627

Incompleteness of available data, caused by end-of-study (administrative) censoring, means that longer term observations beyond the formal follow-up period of the trial are absent. Thus, empirical estimates of survival and CIFs that provide the waiting time distributions will be similarly time-limited. This creates an obstacle in establishing the MGF of components of the MSM, which appears to restrict application of inversion methods. For example, in the LIPID trial, a randomised clinical trial with 5- to 7-year follow-up, there was a high probability (75–80%) of event-free history in the study population during follow-up. Hence, the Kaplan–Meier survival estimator does not descend to 0 at any time during follow-up. Our MSM specifies that patients will experience either stroke or death as first event following randomisation (occurrence of other intermediate events not being considered). Because eventual progression to death is certain, these are complementary first events. However, empirical estimates of CIFs of time from randomisation to these alternative first events (stroke and death) do not sum to 1, even 6 years or more after randomisation. A consequence of the censoring is that empirical estimates of survival distributions and CIFs are incomplete, being undefined beyond the maximum time of follow-up for exit from each state. Hence, nonparametric estimates of the MGFs of the holding times in states are unavailable. Furthermore, the proportions p of patients who will eventually experience stroke before death is unknown. Hence, transmittance through the network is not available without further assumptions. Thus, nonparametric estimation of survival or hazard functions by saddlepoint inversion of transmittance cannot proceed. End-of-study censoring is common in many biostatistics applications (cohorts, trials, etc). Butler and Bronson [3] proposed a nonparametric estimator of transmittances for censoring mechanisms, which are tail-estimable. Their Theorem 3 proves convergence of survival, density and hazard functions provided by saddlepoint inversion of sample estimates to population counterparts. This important general result requires the tail-estimability condition that ‘all exit-time distributions fFij .t /g have intervals of support that are subsets of .0; H /, the support of [censoring variable] C , with H 6 1’. The theorem establishes that a method (rescaling) of redistributing unallocated mass is correct in large samples. However, this condition of tail-estimability is often not met in trials with limited periods of follow-up. Because of this, effects of rescaling are not asymptotically negligible. Rescaling substitutes an estimate of transmittance of the marginal distribution of time to exit from any given state on the basis of the empirical distribution of time to next event [8] . Let O be the largest observed time duration before exit from this state. The rescaled transmittance is that of the conditional distribution given an event in the interval Œ0; . O This approach biases downward the discounted waiting time (the MGF) by omitting holding times that exceed O . Furthermore, without knowledge of the cause-specific hazard functions beyond time O , the proportion of persons eventually following each path is unknown, so the probability of the path and its transmittance are inestimable; see [3, Theorem 4]. Although study subjects will eventually exit state 1 in LIPID with probability 1, many holding times in this state will exceed their duration of study follow-up. Although reallocation of remaining departures to one or other pathway in proportion to numbers observed so far may be plausible, calculating transmittance by substituting for the MGF of the marginal distribution of time to first event, the MGF of a conditional distribution introduces bias. Without application-specific justification, rescaling using transmittances of distributions with support within Œ0; H , with H < 1, provides survival estimates

H. M. HUDSON ET AL.

0.25 0.20

Kaplan−Meier KM conf.limits SP rescaling SP Markov

0.15

Cum.Incidence

0.15

0.00

0.05

0.10 0.00

0.05

Cum.Incidence

0.20

Kaplan−Meier KM conf.limits SP rescaling SP Markov

0.10

0.25

on the basis of biased estimates of the MGF, as the substituted survival time distribution is stochastically smaller than the marginal survival. Any saddlepoint approach to survival estimation then depends on assumptions that are untestable. The problem is best appreciated in the context of simple survival analysis, albeit one in which saddlepoint inversion has few advantages over simpler empirical estimators. Even in the network with no intermediate events, rescaling needs correction, as described in Appendix A. There are many possible ways to deal with this issue of the unallocated mass. One possible solution is to avoid the problem altogether by using fully parametric estimation, but this strategy lacks the flexibility of nonparametric techniques. Another approach is to consider nonparametric estimation in a period of follow-up Œ0; and impose assumptions, which permit suitable extrapolation of survival and CIF distribution estimates thereafter. The quality of the extrapolation can only be judged on its plausibility on the basis of the data at hand and on its effect on estimates within the period Œ0; , and its consistency with data in this time-period. In general, the lack of estimability above the maximum observed survival time makes it impossible to provide a fully nonparametric estimate of the transmittance of time to endpoint that is required by saddlepoint inversion; see [3, Theorem 4]. Figure 4 provided three estimators of cumulative incidence of death (any cause) in LIPID; none of which uses intermediate event data. That is, the network considered has only the single path from randomisation to death (1 ! 3). The time scale origin is the time of randomisation. The figure serves to display the level of accuracy possible in estimating survival functions with the two methods employing saddlepoint inversion as well as confirming the accuracy of the adjustment to scaling given in Appendix A. In this single transition survival context with censoring, saddlepoint estimates use the same information and are directly comparable with the empirical CIF 1 SO .t / derived from the Kaplan–Meier estimator. This is the reference survival estimator shown in Figure 4. Its confidence limits are also shown as dashed lines. Figure 4 also provided two saddlepoint estimators requiring additional assumptions. The first estimator applies a rescaling to estimate transmittance along the sole network path 1 ! 3, before saddlepoint inversion followed by correction. It is described in detail in Appendix A. This appendix includes the definition of the adjusted rescaling we used in Figure 4. The saddlepoint estimator labelled ‘SP Markov’ uses the method developed in this paper (Section 3.2), extrapolating survival beyond 6 years by employing a translated exponential tail (with constant residual hazard). In semiparametric extrapolation, explicit

0

1

2

3

4

5

time on study (yrs)

6

7

0

1

2

3

4

5

6

7

time on study (yrs)

1628

Figure 4. Estimates of cumulative incidence of trial outcome, death of any cause. Solid lines show empirical and two saddlepoint survival estimates by treatment group (left panel: control; right panel: statin). Estimators are (i) the Kaplan–Meier estimator and saddlepoint estimates using (ii) corrected rescaling (red) and (iii) constant hazard extrapolation (blue). Dashed lines show Kaplan–Meier confidence limits.


Statist. Med. 2014, 33 1621–1645

H. M. HUDSON ET AL.

expressions for the CGF and its derivatives (Appendix C) provide the distribution of time from randomisation to the primary endpoint in the network by their use in the saddlepoint density formula. The residual survival extrapolation assumes constant hazard from 6 years onward (death rates 33.8/1000 per year and 24.9/1000 per year in controls and in the pravastatin group, respectively). Provided the correction to scaling is made in this simple survival analysis context, all three estimates of survival up to 6 years shown in Figure 4 are very similar, with estimates based on saddlepoint inversion remaining within the Kaplan–Meier confidence intervals at most times. The correction applied in Appendix A to provide an estimate of survival up to time does not provide an estimate of transmittance (because the MGF remains indeterminate) and so does not generalise to multistate networks. In any network other than a single step to endpoint, we know of no way other than extrapolation to estimate marginal transmittances of component stages. 2.6. Other networks As noted at the conclusion of Section 2.2, explicit expressions for the transmittance (and its derivatives) through any finite-state semi-Markov network are readily calculated from component transmittances. This transmittance is the MGF of time through the system in the case of a single absorbing state (endpoint), as the probability of reaching this endpoint is 1. Thus, the CGF and its derivatives are available in networks with a single absorbing state. Saddlepoint inversion then proceeds using the algebraic identities of Butler and Bronson; see equations (13.12)–(13.14) of [6]. Thus, in cases such as the LIPID randomisation-stroke-death model, it remains only to estimate component MGFs and their derivatives; in Section 3, we describe computations using empirical Kaplan–Meier survival and CIFs, extrapolating when required. We show that the MGFs for competing pathways are readily computable after extrapolation of CIFs, as in equation (10), or the particular case equation (9). For transitions without alternative branches, such as exits from states 2 to 3, extrapolating the Kaplan–Meier survival estimator will provide the corresponding MGF. Alternatively, in any MSM stage, fully parametric distributions may be fitted (corresponding to extrapolation from time D 0). In states with the possibility of re-entry after exit (loops), exit to the state of departure is considered as one of the competing risk events. Saddlepoint inversion is also available with competing risk outcomes (two or more absorbing states). Extrapolation provides proper distributions of times spent before exit from any state, so that the probability of each outcome being realised can be calculated, along with the transmittance to that outcome. We calculated the conditional MGF of Y given that Y < 1, where Y denotes time to a specified endpoint, by dividing the transmittance to that endpoint by the probability of this endpoint being reached, so that the CGF and its derivatives are available. Then we can obtain saddlepoint approximations to the survival and hazard functions, conditional on the specified endpoint, by the methods described in [5]; see Theorem 2.

3. Empirical transform estimates As discussed earlier, in segments of the network, the methods of Butler and Bronson require estimates of the transforms (MGFs) of corresponding survival distributions. In this section, we discuss empirical estimation of these component MGFs in segments where the transition distribution is to a single state, or alternatively, transition can occur to one of several competing successor states. Observed survival times then need to be accompanied by cause of the exit from the current state (or indicator of censoring). 3.1. Kaplan–Meier integrals and truncation


Statist. Med. 2014, 33 1621–1645

1629

Let X be the waiting time to the first observed event and H be a fixed upper bound on followup (censoring times) in a network of progressive states. We seek empirical estimators of the MGF M.s/ D E exp.s X / of time to first event. PIn the absence of censoring, with a random sample of n observations, the empirical estimator exp.sXi /=n then provides the necessary component MGF. With censored observations present, a corresponding estimator is the Kaplan–Meier integral [14, 15] of s .x/ D exp.s x/, given as Mn .s/ D R1 0 exp.s x/ dSn .x/, where the empirical survival distribution Sn is the Kaplan–Meier estimator. The empirical Kaplan–Meier integral Mn .s/ is strongly consistent, and converges ‘in mean’, as n ! 1, to its population counterpart. Conditions on censoring for convergence to this limit can be found in Reference [15]. These conditions preclude tail-dominating censoring schemes (like type I censoring).

H. M. HUDSON ET AL.

However, when the censoring time distribution has bounded support Œ0; H , convergence is to a different limit, a truncated Kaplan–Meier integral. For any 6 H , the right extreme of the distribution of censored survival time, write X D X 1fX6 g C . C R/ 1fX> g D X C R 1fX> g

(1)

where 1A is the indicator of occurrence of an event A, X D X 1fX6 g C 1fX> g

(2)

and R D X is the residual survival after time . The random variable X is the survival time, truncated at . By Stute’s Kaplan–Meier integral theory [15, equation (2.7)], its MGF is consistently estimable. Moreover, its MGF, M .s/ D E Œexp.s X /, is related to that of X , by E Œexp.s X / D E Œexp.s X / S. /e s f1 E Œexp.s R/ j R > 0g Any censoring prior to in X is allowed for in calculating the first term of the right hand side using its consistent Kaplan–Meier integral estimate. This Kaplan–Meier integral of survival truncated at exists for all s, because of the bounded limits of the integral. It follows from the aforementioned expression that the MGF of X can be recovered from this first term by subtracting the second term, involving the MGF of the distribution of residual survival, beyond . Thus, when knowledge of the empirical survival distribution is restricted to the domain Œ0; , the MGF is not estimable without additional assumptions. By specifying a parametric distribution of residual survival among survivors at follow-up time for the second term, corrections for limited trial follow-up can be implemented, using the aforementioned expression. Such extrapolation beyond a chosen point 6 H is necessary when time to the first event has non-negligible probability that X > H . Suzukawa [16] extended Stute’s proof of strong consistency to CIFs of competing risks, and [3] extended Suzakawa’s work to SMP network states. In SMPs, the number of individuals entering the state is a random variable. In both cases, an upper limit H , in SMPs specific to each newly entered state, is important to the convergence. Extrapolating residual survival in each successive network stage following first event (with time-scale beginning from entry to each new state and extrapolation beyond 6 H ) will then provide estimates of the transmittance in each stage. 3.2. Extrapolating survival with bounded period of follow-up Suppose a reliable estimator of the survival function to first event, S.t / D P .X > t /, is available for follow-up times 0 6 t 6 . Define residual survival beyond as random variable R D X . Working from the estimate of S. /, we may then approximate survival in the domain t > by extrapolating the function S.t / beyond time . A common extrapolation of the Kaplan–Meier estimator employs translated exponential tails. The extrapolation of S.t / is then defined as S.t / D S. / exp Œ .t / for t > . Conditional on R > 0, residual survival R D X beyond is governed by the negative exponential distribution with constant hazard for which the cumulative hazard function is ƒ0 .u/ D u for u > 0, so that S.t / D S. / expŒƒ0 .t / for t > . The linear hazard rate distribution (LHRD) [17, 18], also known as the linear failure rate distribution, is a family of survival distributions including both negative exponential and Rayleigh distributions as special cases. In empirical analysis for an earlier paper, we found that the LHRD provided a better fit to observed survival of participants in the LIPID trial to first event [1]. The LHRD assumption for residual survival beyond trial time horizon specifies S.t / D S. / exp ˛ .t / .t /2 =.2ˇ 2 / D S. / expŒƒ0 .t /

1630

for t > , where now the cumulative hazard ƒ0 .u/ D ˛u C u2 =.2ˇ 2 / with ˛ > 0; ˇ > 0, for u > 0, is quadratic. Copyright © 2013 John Wiley & Sons, Ltd.

Statist. Med. 2014, 33 1621–1645

H. M. HUDSON ET AL.

More generally, extrapolation extends S.t /, defined on t < , by choice of any non-negative, increasing function ƒ0 .:/ on Œ0; 1/, which increases without bound. Equivalently, specify the survival function S0 .t / D exp Œƒ0 .t /. The extrapolation is defined by S.t / D S. / exp Œƒ0 .t / D S. / S0 .t /, for t > , so that the conditional distribution of R D X , given R > 0, then has survival function S0 .r/ and cumulative hazard function ƒ0 .r/ for r > 0. 3.3. Transitions with competing risks Consider two competing risks, and empirical estimation of the CIF. Events are assigned two exclusive causes—1 and 2—depending on successor state in the network, define the CIF of cause 1 as follows: Z t P .X > u/ 1 .u/ du (3) F1 .t / D 0

where X is the waiting time to the first event and 1 .t / is the cause-specific hazard function for events with cause 1 [19, Section 8.2.3]. For example, in the LIPID randomisation-stroke-death model, exit from state 1 will be either due to death (cause 1) in the transition from state 1 to 3 or stroke (cause 2) in the transition from state 1 to 2. Empirical estimation of each CIF employs the Kaplan–Meier estimator of the survival function of X and the empirical proportions of events at risk points, which are of the specified cause. Thus, for distinct ordered times of events t1 < : : : < tm (not including censored times), FO1 .t / D

X i jti . Extrapolation now has two components: (i) extrapolation of the survival function, as described in Section 3.2; and (ii) specification of the hazard rate of events of each cause beyond time . Now for cause 1, from equation (3), Z t F1 .t / D F1 . / C S.u/ 1 .u/ .u/ du 0

(4)


Statist. Med. 2014, 33 1621–1645

1631

for t > . Because, by definition, .t / D S .t /=S.t /, we have Z t S 0 .u/ 1 .u/ du F1 .t / D F1 . /

H. M. HUDSON ET AL.

Thus, for t > , dF1 .t / D 1 .t / dS.t /

(5)

3.6. Proportional incidence assumption A simple form of extrapolation will utilise constant relative hazard of cause 1 in equation (5), which we term a proportional incidence assumption. Let .t / be the hazard function of time to first event of any cause. Suppose the proportion of events of type 1, 1 .t /=.t /, is independent of time, beyond time . Then 1 .t / D , where is constant, for t > . Proportional incidence of events of different causes will apply after time . Then, by equation (4), F1 .t / D F1 . / C ŒS. / S.t / D F1 . / C S. / f1 exp Œƒ0 .t /g

(6) (7)

for t > and any specified extrapolation model, that is, specified non-negative increasing function ƒ0 .u/ on u > 0 increasing without bound. This expression simplifies extrapolating CIFs of competing events with proportional incidence. In the LIPID randomisation-stroke-death model, with the constant proportion of stroke among first events in subjects not experiencing a first event before time , the probability p of stroke among all subjects is p D limt !1 F1 .t / D F1 . / C S. /; from equation (6). Thus, the proportion of subjects p who experience stroke as first event may be estimated by an adjustment to the empirical estimator of the CIF for this cause at time . When is F1 . /=ŒF1 . / C F2 . /, as in rescaling, it is easily shown that p D . Equation (4) provides a general method of extrapolating CIFs, not dependent on proportional incidence assumptions. For example, when relative hazard rate is not constant, but the hazard of first event is constant after time , we can calculate CIFs using Laplace transforms—for events of cause 1, the Laplace transform of 1 .t /—as can be seen from this equation. 3.7. Moment-generating function of waiting time to first event, given its cause The empirical CIFs of waiting times to events of each cause also provide empirical transmittances as non-parametric estimators. The transmittance for the sub-distribution (CIF) of cause 1 is defined by Z

1

T1 .s/ D

exp.s t / dF1 .t /

(8)

0

where F1 is the sub-distribution of time to event of cause 1. Note that defining transmittance for first event following randomisation in this manner in the randomisation-stroke-death network provides T1 .s/ D T12 .s/ D p M12 .s/ where M12 is the MGF of conditional distribution of time to stroke (cause 1) as first event, given this occurrence. Similarly, T2 .s/ D T13 .s/ D .1 p/ M13 .s/ is the transmittance for cause 2 (death as first event). Note that the probability p is implicit in the definition of the aforementioned transmittance, once the sub-distribution is specified. It is also important to note that the estimator of transmittance employing the empirical CIF in place of F1 in equation (8) is undefined in the case of type 1 censoring as the CIF is indeterminate for times beyond . Instead, we propose the use of a semiparametric estimator, which substitutes the transmittance of a composite of the empirical CIF and its parametric extrapolation beyond a selected time horizon . From equation (5), the transmittance T .s/ of the improper distribution (CIF) of cause 1 is Z

Z

1

0

Z

exp.s t / dF1 .t / D

1

exp.s t / dF1 .t / 0

exp.s t / 1 .t / dS.t /

(9)

which we write T .s/ D T o .s/ C T e .s/, for convenience. Now, with proportional incidence beyond , and extrapolation (Section 3.2) of first-event survival S.t / D S. / S0 .t / for t > ,

1632

e

Z

1 s

T .s/ D

Z

1

exp.s t / 1 .t / dS.t / D e S. /


exp.s u/ dS0 .u/

(10)

0

Statist. Med. 2014, 33 1621–1645

H. M. HUDSON ET AL.

Empirical estimates of the CIFs are obtained by the method of Kalbfleisch and Prentice and reviewed in [2]. Implementations are readily available in software, for example, in packages survival and cmprsk in R. With F1 .t / estimated by the empirical CIF for cause 1, FO1 .t /, a step function, the first integral T o .s/ on the right hand side of equation (9) reduces to a finite sum, readily computable. In equation (10), an extrapolation of survival past is needed. As appropriate choice of distribution, we recommend a simple parametric form, for example, using exponential or LHRD for residual survival S0 .u/. The Appendix contains the MGFs of the residual survival distributions for these extrapolations.

4. Application in randomised controlled trials There have been few applications of MSMs reported in statistical or clinical medicine journals since the pioneering evaluation of Andersen et al. [21] using three-state models. Methods for analysis of a (uncontrolled) cohort study employing a multi-stage model of states in burns patients were described in [22], whereas clinical reporting of recurrent event data with terminal event death using an MSM appears in an analysis of metastatic cancer [23]. In previous work [1], we described application of saddlepoint inversion using parametric survival distributions in network components, and accounting for censoring. The method was applied to the LIPID study data in the simple illness-death model informed by the intermediate event of initial stroke and terminal outcome of all causes of death. Comparison was made with standard survival approaches and HR findings. In LIPID, we found the time to first event among competing events, initial stroke and death of any cause, to be poorly fitted by a negative exponential distribution, but well fitted by the LHRD. Transform inversion of time to the study endpoint showed excellent agreement with Kaplan–Meier estimator ignoring the intermediate event of stroke, and indicated a previously unreported cumulative benefit of pravastatin. In the presence of this evidence of non-PHs, functionals for randomisation-based inference demonstrated statistically significant benefits of statin therapy. In the following section, we conduct survival and HR estimation by saddlepoint inversion of firstpassage transmittance estimated using our semiparametric extrapolations within the randomisationstroke-death network of LIPID. This approach reduces the need to correctly specify parametric components of MSMs when estimating treatment effects in trials with intermediate events affecting outcomes. We calculate first-passage transmittance from transmittances of individual components of the network. In order to define transmittance within each component, we calculated the transmittances associated with the exit from a current state using empirical distributions [7, 8, 24]. We extended empirical distributions by parametric distributions of residual survival after a specified period of follow-up. The study served to demonstrate the feasibility of integrating MSM elements within a semi-Markov network. In the LIPID study, stroke was an intermediate outcome in an illness-death SMP model of survival. Most study participants remained alive and stroke-free at study end (with pathway to final study endpoint remaining unknown). In this section, we study the suitability of transform methods in estimating and comparing survival functions in the treatment (pravastatin) and control (placebo) arms of the LIPID trial. 4.1. Results: survival to first event


Statist. Med. 2014, 33 1621–1645

1633

We showed empirical survival and CIF distributions of time to first event (stroke or death) following randomisation in the LIPID study in Figures 5 and 6. Note that substantial end-of-study censoring occurs after 5 years and before 7 years following randomisation (Table II). Listing proportions lost to follow-up by years of follow-up, we chose the cutpoint to be 6 years (denoted 1 in Figure 2) following randomisation as the time-horizon appropriate for empirical estimation of survival to first event. After this time, empirical estimates of survival and CIFs were considered less reliable because of low event numbers and end of follow-up at trial closure. We based LHRD parameter estimates used in the extrapolation component M e .s/ of the MGFs M12 .s/ and M13 .s/ on times to first event, or censoring, of all study subjects still at risk 2 years after study entry, unless otherwise noted. LHRD parameters for time to first event were fitted by maximum likelihood.

0.20

H. M. HUDSON ET AL.

0.05 0.02

Cum.Incidence

0.10

placebo pravastatin

1

2

3

4

5

6

time on study (yrs)

Figure 5. Cumulative incidence of first events. The log scale shows greater reductions in hazard ratio of first event among pravastatin-assigned patients relative to controls with increasing time on study.

0.12

0.14

Cumulative event incidence to T=6 years

0.08 0.06 0.00

0.02

0.04

Cum.Incidence

0.10

death − placebo stroke − placebo death − pravastatin stroke − pravastatin

0

1

2

3

4

5

6

time on study (yrs)

Figure 6. Cumulative first-event incidence by treatment group and type of event (stroke, death). Table II. Numbers at risk of first events by randomised treatment arm and years since study entry. Years since entry

Control Statin

0

1

2

3

4

5

6

4502 4512

4411 4422

4292 4328

4189 4237

4041 4127

3660 3750

1706 1809

1634

To apply the proportional incidence assumption, we specified a common parameter, D 0:76, of proportion of deaths among strokes and deaths in both pravastatin and control arms. We based this choice on observed proportions in data on all study subjects, and those with survival exceeding 2 and 5 years (Table III). Copyright © 2013 John Wiley & Sons, Ltd.

Statist. Med. 2014, 33 1621–1645

H. M. HUDSON ET AL.

Table III. Observed proportions of deaths among first events. Year of first event By time interval

0–1

1–2

2–3

3–4

4–5

5+

0:77 0:74

0:71 0:76

0:73 0:70

0:76 0:79

0:77 0:82

0:81 0:72

0+

2+

5+

0:761 0:758

0:770 0:761

0:806 0:723

Control Statin After specified follow-up Control Statin

Time through the system in the LIPID trial also depends on the survival distribution S23 of patients who experienced a stroke as their first event. The intermediate stroke state here included all initial strokes following randomisation. We assessed time to stroke for any first stroke, whether fatal or not. In a previous parametric analysis we conducted, it had been necessary to redefine the stroke state to include only non-fatal strokes (defined as those in whom death did not occur within the subsequent 14 days) as intermediate events, so that a simple constant hazard model of subsequent death would provide a satisfactory fit to the data. Such redefinition of states is not entirely satisfactory, and not necessary in semiparametric modelling. To estimate the survival distribution S23 , we use the Kaplan–Meier estimator to 24 months following initial stroke, extended by assuming constant hazard of death thereafter, an exponential distribution of residual survival. (Constant hazard is a limiting case of linear hazard rate assumptions. A more general LHRD of residual survival did not appear to be warranted here). This approach captures high early mortality following stroke (during the 24-month period) and subsequent low event rates following recovered stroke. 4.2. Integration: time through the system.

0.95 0.90

KM−placebo KM−pravastatin SP−placebo SP−pravastatin P−placebo P−pravastatin

0.85

Survival probability

1.00

Applying an LHRD of residual survival to first event, and proportional incidence assumption beyond 6 years, together with the survival from stroke to death estimated as earlier, we obtain the survival and, hazard to final study endpoint, any cause of death. The fitted survival function is presented with label ‘SP’ in Figure 7, where comparison is made with the empirical survival estimator (Kaplan–Meier estimator, ‘KM’) and with a fully parametric saddlepoint estimator (‘P’). We presented corresponding estimators of cumulative incidence in Figure 8.

0

1

2

3

4

5

6

1635

time on study (yrs)

Figure 7. Survival estimates: empirical, parametric and semiparametric.


Statist. Med. 2014, 33 1621–1645

0.10

0.15

KM−placebo KM−statin SP−placebo SP−statin

0.00

0.05

Cum.Incidence

0.20

0.25

H. M. HUDSON ET AL.

0

1

2

3

4

5

6

7

time on study (yrs)

Figure 8. Fitted cumulative incidence using linear hazard rate distribution extrapolation to residual survival beyond 6 years. Linear hazard rate distribution parameters are estimated among those remaining at risk at 2 years.

0.15

Probability

0.20

0.25

0.30

death stroke extrap: M1 extrap: M2 LHRD

0.00

0.05

0.10

0.15 0.00

0.05

0.10

Probability

0.20

0.25

0.30

death stroke extrap: M1 extrap: M2 LHRD

0

5

10

Years

15

0

5

10

15

Years

1636

Figure 9. Empirical cumulative incidence functions of first-event outcomes all-cause death and stroke with parametric extrapolations from 6 to 15 years in control (left panel) and pravastatin (right panel) treatment groups. Parametric extrapolation assumes: (a) M1: constant hazard (Markov) extrapolation from T=6 with event rate and treatment specific proportional incidence estimated in those event free at 5 years; (b) M2: extrapolation as for M1 but with event rate estimated beyond 4 years and constant relative hazard of death 0.76; (c) linear hazard rate distribution extrapolation with parameters estimated among those event-free at 4 years and constant relative hazard of death 0.76.


Statist. Med. 2014, 33 1621–1645

H. M. HUDSON ET AL.

We have shown the proportional incidence model with LHRD residual extrapolation of Kaplan–Meier survival to first event in Figure 9, and observed proportions of deaths among first events appear in Table IV. 4.3. Integration: statistical measures and inference We estimated the hazard function for the survival time in each study arm by saddlepoint inversion of the CGF of time through the network, using empirical estimates of the network components (transmittances). The HR is represented (as bold curve) as a function of time since randomisation in Figure 10. By repeating the saddlepoint inversion, in each data sample re-estimating the parametric extrapolation among those remaining at risk at 4 years, we calculated 2500 bootstrap replicates of the HR at monthly intervals. Using the bootstrap percentile method, the 2.5% and 97.5% quantiles of bootstrap estimates provided a 95% confidence range for the HR at each time shown as a shaded area within this Figure. For purposes of comparison, we have shown the fully parametric fit and bootstrap confidence intervals in the right panel of Figure 10. The trial reported HR 0.78 for treatment allocation to statin, shown as a horizontal line on each panel of the figure together with the associated 95% confidence range.

Table IV. Observed proportions (prop) of deaths among first events. Controls

Pravastatin

Year

Deaths

Strokes

Prop

Deaths

Strokes

Prop

0–4 4–5 5–6 6+ All

341 103 105 40 589

119 31 28 7 185

0:74 0:77 0:79 0:85 0:76

289 98 63 23 473

96 22 22 10 150

0:75 0:82 0:74 0:70 0:76

Parametric 1.2 1.0 0.8

Hazard Ratio

0.4

0.6

0.8 0.4

0.6

Hazard Ratio

1.0

1.2

Empirical

1

2

3

4

5

6

7

1

2

3

4

5

6

7

years on treatment

Figure 10. Time dependent hazard ratio of treatments.


Statist. Med. 2014, 33 1621–1645

1637

years on treatment

H. M. HUDSON ET AL.

1.5 1.0 −0.5

0.0

0.5

Cumulative Survival Gain (months)

1.0 0.5 −0.5

0.0

Cumulative Survival Gain (months)

1.5

2.0

Parametric

2.0

Empirical

1

2

3

4

5

6

7

years of follow−up shaded: 95% confidence

1

2

3

4

5

6

7

years of follow−up shaded: 95% confidence

Figure 11. Net treatment benefit from statin therapy.

1638

Because HR is not constant over time, we estimated net treatment benefit, a functional of time T estimating expected survival gain on pravastatin treatment from years 1 to T, a measure proposed by Lagakos [25]. We provided estimated net treatment benefit, the difference in cumulative incidence of any cause of death, in Figure 11, together with (shaded) confidence ranges derived from bootstrap replicates of this statistic. We computed permutation tests of the statistical significance of treatment differences: (i) using differences in survival estimated (by saddlepoint) at 5 years; and (ii) using net survival benefit from 1 to 5 years or from 1 to 6.5 years. we calculated approximate p-values in (i) and (ii) using the proportion of permutation resamples for which the statistics (differences between assigned groups in survival at 5 years or net survival benefits) estimated as great a differences between treatment groups as that observed in the study. The numbers of resamples with differences as large as that observed as a fraction of the total number and the approximate p-values were (i) 1/2500 (permutation test p < 0:001); and (ii) 1/2500, 0/2500 (each with permutation test p < 0:001). We calculated a bootstrap test of non-PHs using 10 000 bootstrap resamples, in each data sample re-estimating the parametric extrapolation among those remaining at risk at 4 years and calculating differences in log HR estimated at 1 and 5 years. We used the bootstrap percentile method to obtain the sampling distribution of the difference in log HR between its values at 1 and 5 years. An approximate p-value for the test of non-constant HR is then twice the area under this distribution below (difference) 0. Corresponding numbers of bootstrap resamples and approximate p-value were 2 218=10000.P D 0:04/. We achieved higher significance of net survival benefit by extending the truncation time of benefit to 6 or 7 years. Permutation tests of treatment differences using differences in survival at 5 years and from use of net survival benefit are therefore comparable with the all-cause mortality logrank p < 0:001. Neither tests of treatment differences in survival at specified times nor tests of net survival benefit of treatment depend on the PHs assumption, which the evidence from comparison of log HR at 1 and 5 years rejects. This finding of accumulating benefit with time since randomisation is sharper than would be indicated by the non-significance of commonly used tests for this interaction, which make no use of intermediate event data in survival analysis. The 95% bootstrap confidence interval for difference in Copyright © 2013 John Wiley & Sons, Ltd.

Statist. Med. 2014, 33 1621–1645

H. M. HUDSON ET AL.

log HR between 1 and 5 years was Œ0:519; 0:008 with mean 0:260 corresponding to an expected 23% increase in pravastatin’s effect on mortality reduction after 5 years compared with that after 1 year on treatment. 4.4. Diagnostics and model sensitivity

0.25

0.25

We can provide an empirical estimator of survival not dependent on transform methods. This estimator derives from the adaption of the Aalen–Johansen estimator [26] to SMPs; see Example X.1.7 of [27]. The comparison with saddlepoint inversion of transforms of empirical component distributions provides a diagnostic examining the accuracy of saddlepoint approximation in this model. We obtained the Aalen–Johansen NPML estimator for outcome all-cause mortality in the simple randomisationstroke-death LIPID model using the R package mstate [28]. The resulting estimates of cumulative incidence were very close to the Kaplan–Meier estimates shown in Figure 8. Therefore, it appears that allowance for intermediate events data has little effect on the survival distributions estimated (for pravastatin and control groups) using the Aalen–Johansen estimator. Although benefit of the use of intermediate event data was mainly the improved significance of tests of treatment comparisons and interaction (non-PHs), this finding allows the use of the Kaplan–Meier estimator as a proxy standard against which the accuracy of saddlepoint estimators can be assessed (see results in the succeeding text). Using Table IV, we can evaluate the constant relative hazard assumption employed in saddlepointbased inference of LIPID data: (i) for differences in relative hazard between treatment and control; and (ii) for differences over time. A chi-square test of differences between observed and fitted counts in the table provides 2 D 3:71 for controls, 2 D 3:18, each non-significant on 3 df. We obtained similar results for tests of trend (with time interval) in the table. Thus, the data are consistent with the assumption of constant relative hazard.

Cum.Hazard

0.15

0.20

Kaplan−Meier KM conf.limits SP extrap M2 SP extrap LHRD

0.00

0.05

0.10

0.15 0.10 0.00

0.05

Cum.Hazard

0.20

Kaplan−Meier KM conf.limits SP extrap M2 SP extrap LHRD

0

1

2

3

4

5

time on study (yrs)

6

7

0

1

2

3

4

5

6

7

time on study (yrs)


Statist. Med. 2014, 33 1621–1645

1639

Figure 12. Fitted cumulative hazards: time to Long-Term Intervention with Pravastatin in Ischaemic Disease endpoint (all-cause death). Kaplan–Meier estimates and saddlepoint estimates based on constant hazard extrapolation (M2) and linear hazard rate distribution extrapolation by treatment group (left panel: control; right panel: statin.)

H. M. HUDSON ET AL.

We illustrated sensitivity of the estimated survival function to specification of the extrapolation by fitting constant hazard residual extrapolations in place of the linear hazard rate residual extrapolations of cumulative incidence of first event. Extrapolations were from 6 years; models employed constant relative hazard 0.76 of death among first events, common to pravastatin and control groups, except in the case of the Markov extrapolation model M1. In model M1, more optimistic assumptions of reduced events among patients randomised to receive pravastatin—suggested by the observed event counts of subjects remaining event-free past 5 years—specified relative hazard 0.723 of death among first events in the pravastatin group, but 0.806 in the control group. Markov extrapolation models M1 and M2 specified constant hazard of first event; rates were estimated for statin and control treatment groups from event rates among patients event-free beyond 4 years. Note that a constant hazard model is a poor fit to first-event data for the interval from 0 to 6 years, where an LHRD provides a statistically significant improvement. Nevertheless, with the paucity of exposure beyond 6 years, and the apparent increase in event rates between 5 and 6 years, the assumptions of extrapolation models M1 and M2 represent plausible scenarios for later years’ experience of first event. Estimates using extrapolation model M1 are similar to those using model M2, shown in Figure 12. To simplify viewing, we have not shown estimates for model M1 are not shown in Figure 12. Figure 12 indicates that the saddlepoint survival estimates are in close agreement for a variety of extrapolation assumptions. Saddlepoint estimates using LHRD residual extrapolation are in close agreement with Kaplan–Meier estimates, a proxy for the NPML estimates, but are smooth. It may be observed in the figure that an assumption of constant hazard in first event beyond 6 years reduces the accuracy of the saddlepoint solution in earlier years, which estimates a cumulative survival hazard above that of the empirical estimator (the negative logarithm of the Kaplan–Meier estimate of survival probability), particularly between 2 and 6 years after randomisation. However, these estimates based on constant hazard extrapolation remain within the Kaplan–Meier confidence bands, skirting the upper limit of that attributable by chance as indicated by this model-free estimation of cumulative hazard.

5. Conclusions

1640

In trials, a semi-Markov model more fully utilises information from a network of intermediate events than does a standard survival analysis. The nonparametric approach of [3] provides an effective method of estimating transmittances for successive transitions between states of the model using intermediate event data. Simple algebraic calculations allow the net transmittance along all paths from entry to the final endpoint to be estimated. Saddlepoint inversion of this transmittance provides survival and hazard estimates of survival time to endpoint. This approach removes the need to specify suitable parametric representations; empirical transforms are substituted instead in saddlepoint calculations. However, the method requires the transmittance to be rescaled should the longest follow-up time to next event be censored. Its properties depend on the assumption that survival to next event be tail-estimable, given censoring. The tail-estimability condition is not met in most clinical trial applications where limited trial follow-up causes the truncation of Laplace–Stieljes integrals necessary for saddlepoint inversion. In such situations, we propose a semiparametric estimation of the transmittance, replacing the rescaling but requiring extrapolation of residual survival and proportional incidence assumptions for transitions with competing risks. Our approach adds to the flexibility of saddlepoint approximation in MSM networks. We have shown that allowance for a truncated period of trial follow-up by extrapolation can reduce the bias inherent in saddlepoint approximations employing rescaled transmittances. Residual extrapolation provides a transmittance of closed form, which serves as a ready alternative to transmittance on the basis of individually specified parametric models. In controlled clinical trials, the method becomes attractive as a simple yet flexible semiparametric method of studying treatment effects and interpreting net effects on survival and HRs in MSMs. Our approach explicitly specifies extrapolation assumptions, which may affect survival comparisons. Different assumptions about residual survival may then be considered in the MSM and sensitivity analysis of findings evaluated. Extrapolation of survival and hazard functions from randomisation to a selected trial endpoint are also available with this approach. The use of saddlepoint inversion makes valid inferences using bootstrap inversion feasible and practical. These methods supplement MSM analysis in randomised trial applications in a synthesis estimating survival functions and HRs by fully utilising information gathered on intermediate events that precede the trial terminal outcome(s). This analysis is not limited by network complexity. Copyright © 2013 John Wiley & Sons, Ltd.

Statist. Med. 2014, 33 1621–1645

H. M. HUDSON ET AL.


Statist. Med. 2014, 33 1621–1645

1641

This paper also provides an interpretation of MSMs from a clinical perspective through integration of MSM components in making findings for survival to primary outcome. The LIPID study, a randomised trial important in establishing the beneficial effects of statin therapy in a large sample of known CHD patients at high risk of CHD events such as cardiovascular death, acute MI and stroke, has been used as illustration. The initial study showed that pravastatin therapy over 6 years reduced CHD and all-cause mortality and cardiovascular events in this population, results later confirmed in many other studies including the 8-year LIPID follow-up [29]. A parametric MSM using stroke as intermediate event (i.e. a simple illness-death model) provided evidence of a cumulating benefit of treatment with time, not consistent with a fixed proportional HR model; see [1]. From a clinical viewpoint a PH model may only provide an average treatment effect in the period of follow-up which underestimates the true treatment effect if the intervention under investigation offers greater benefit for extended follow-up as recently observed in a meta-analysis of statin therapy [30]. The parametric approach in [1] provided excellent fits to the empirical survival to the endpoint (all cause mortality) in both arms. The present paper demonstrates that, in the LIPID randomisation-stroke-death model, the fit of semiparametric estimators, utilising empirical CIFs in place of parametric distributions of time to next event, achieves the same excellent fit and confirms the cumulating benefit of statin therapy under reduced assumptions. The improved HR estimates and increased statistical significance of tests of treatment effect provide evidence of the information gain from use of intermediate event data in LIPID. Further analyses of the LIPID data using MSMs could examine key pre-specified cardiovascular disease outcomes, including acute myocardial infarction in addition to stroke as intermediate endpoints. Longterm follow-up subsequent to the trial provides the opportunity, in future work, to validate extrapolation of survival functions estimated from our analysis with empirical estimates from the follow-up study of residual survival in the same study groups. Butler and Bronson drew attention to the potential for bias in saddlepoint solutions based on rescaled transmittances of empirical distributions with unallocated mass. We evaluated sensitivity of fitted allcause survival distributions within the LIPID trial to the extrapolation assumptions by varying proportional incidence assumptions by treatment group and fitting a variety of residual survival models. In competing risk and simple illness-death networks, Nelson–Aalen and Aalen–Johansen empirical survival estimators are available that allowed us to check the accuracy of saddlepoint approximations in the period of follow-up of the study. Sensitivity to the extrapolation employed in saddlepoint can therefore be evaluated in simple models. Our preferred extrapolations are those most consistent with available data for the time interval under study, yet parsimonious in parameters. Thus, LHRD extrapolation was preferred to constant hazard extrapolation of first event after 6 years, whereas constant hazard extrapolation of residual survival beyond 2 years after stroke was assumed. Our evaluation of two constant hazard extrapolations of first-event incidence beyond 6 years indicated mild variation in saddlepoint survival estimates, though still within the range of chance. However, when extreme departures in future survival distribution from the observable pattern were specified in modelling time to first event—in analyses not presented here—the lack of smoothness in the transforms that results led to a consequent poor approximation from saddlepoint inversion. This sensitivity means that the choice of extrapolation is important to the saddlepoint methodology. In many cases, it rules out ad-hoc solutions such as rescaling by redistribution of mass, which assigns no probability of events after the follow-up period of the study. Equally, it rules out Efron’s redistribute-to-the-right (RTR) rule in the form in which those remaining event-free past the maximum observed event time are assumed to experience an event shortly after [31]. Although indirect estimation of SMPs via saddlepoint-based transform methods has many advantages, it is not flawless. As pointed out by a referee, the saddlepoint density approximation can be poor in some cases, particularly when bimodal distributions must be fitted [32]. Quantities making use of the density approximation include hazard functions, densities and subdensities; generally, saddlepoint approximations of smoother quantities such as a survival function or CIF are accurate. Circumstances leading to poor approximation are those where the true density has additional high frequency components, providing difficulties for any Fourier approximation method [33]. Fortunately, such circumstances are unlikely to occur in survival analysis, where many distributions are unimodal with densities close to families for which the saddlepoint approximation is excellent (Gamma, normal, inverse Gaussian) so that density estimates providing HRs are likely to be very accurate. Multistate modelling of recurrent event data using saddlepoint methods can be more challenging [11]. Further issues arise when data suggest that relative hazard of some competing risks is not constant. As noted in Section 3.6, extrapolation of CIFs can proceed without the assumption of proportional incidence.

H. M. HUDSON ET AL.

Appendices Appendix A. Survival analysis as a network with one transition Consider a standard survival analysis context with independent right censoring and random variables .Y; ı/. Here, survival time Y is non-negative with distribution G.y/; for y > 0, and failure status ı takes values in f0; 1g. A random sample of n observations from .Y; ı/ provides the sample .ti ; ıi / for i D 1; : : : ; n and associated statistics O D t.n/ , the maximum observed time at risk, and y0 , the maximum observed failure time. Kaplan–Meier estimator extrapolation If individuals remain at risk at O , the Kaplan–Meier estimator is undefined past this maximum observed period of follow-up [34]. A number of commonly applied approaches exist for extrapolation when an estimator of the survival function is required over an extended range, as in the estimation of mean survival time and Kaplan–Meier integrals. Methods providing a proper survival function include as follows: (i) parametric estimation of survival; (ii) truncating survival at the largest observed survival t.n/ (Efron’s method); and (iii) assuming translated exponential tail survival beyond a truncation point. Refer to [35], Section 4.2, Practical Notes 2 and 3, for further details. It is sometimes desirable to further restrict the Kaplan–Meier estimator to a domain on which the estimator is more reliable [36]. For this purpose, we permit choice of an upper limit < O at which point extrapolation commences. Scaling in survival analysis In the case of a network comprising a single stage to endpoint (e.g. the transition 1 ! 2), consider saddlepoint estimation of cumulative incidence G.t / (the complement of the survival function). The transmittance (or MGF, because the path probability is 1) for this single stage is Z

1

T .s/ D

exp.st /dG.t / 0

[3] defined an empirical estimator of transmittance employing rescaling; see Section 4.2, equation (3). It provides O O1 /1 TO .s/ D G.

Z

1

O / exp.st /d G.t 0

where 1 GO is the Kaplan–Meier estimator and O1 D t.n/ is the largest observed survival time (including censored times). O / This definition of rescaling leaves unresolved the definition of the integral beyond t.n/ , because G.t is formally undefined for t > t.n/ [34]. For 6 t.n/ , define O /1 TO .s/ D G.

Z

O / exp.st /d G.t 0

1642

the transmittance of the conditional empirical survival distribution given survival not exceeding . Then O / remains TO , for the choice D t.n/ , accords with the empirical estimator of transmittance cited if G.t constant for t > (an improper survival distribution corresponding to Gill’s extension [37] of the Kaplan–Meier estimator). O / = G. O / is a proper distribution function on Œ0; , so Using Gill’s extension, PO .Y 6 t jY 6 / D G.t O O T .s/ is an MGF. Furthermore, T .s/ when defined this way with D t.n/ satisfies the property ascribed to rescaling by [3], its equivalence with a RTR algorithm. In this algorithm, the exit simulation is performed again in its entirety (a rejection step in the simulation) if t.n/ is a censoring time and the algorithm stops at t.n/ . Copyright © 2013 John Wiley & Sons, Ltd.

Statist. Med. 2014, 33 1621–1645

H. M. HUDSON ET AL.

Redistribute-to-the-right algorithm properties Note that the imposition of a rejection step in Efron’s RTR algorithm implies that the survival random variable Y generated is bounded above by y0 , the maximum observed time to outcome. As a consequence of the rejection step, the variable generated is stochastically smaller than a similarly generated variable produced using RTR but with completed data. In these completed data, those individuals previously censored after y0 continue to be followed to the endpoint under study. Stochastic dominance follows because each RTR simulation on completed data matches that on the original set except for completed censored observations whose survival will exceed t.n/ . But the RTR algorithm with completed data determines the Kaplan–Meier weights providing a consistent estimator of the marginal survival distribution G.t / [31]. As a consequence of stochastic dominance, the MGF of the variable generated using the RTR algorithm with the rejection step is a biased estimator of the MGF of the marginal distribution G. Other approaches might be appropriate in using the redistribution-to-the-right algorithm. For example, observations that lead to rejection could be retained with an additional survival period imputed if an assumption is made specifying the distribution of residual survival. A particular case is Efron’s assumption , where the residual survival is set to 0. It is important to note that rescaling as described in [3] estimates the survival time Y conditional on Y 6 O , whose distribution has support Œ0; y0 (and thus MGF existing everywhere on the real line). As a consequence, saddlepoint inversion of estimated transmittance approximates a conditional distribution of Y . Correcting rescaling Without correction, this scaled estimator approximating the conditional distribution is a poor approximation to the marginal survival time distribution. Specifying cut-off time (e.g. D 6 years), we correct for censoring beyond this cutpoint as follows. (Without such correction, the rescaling method provides a clearly erroneous result, which is not displayed in Figure 3). To estimate the marginal distribution of survival time, elementary probability calculations provide h i O / 1 PO .Y > tjY 6 / PO .Y 6 t / D G. (11) O / is the probability of survival not exceeding years, using the Kaplan–Meier estifor t < , where G. O mator, and P .Y > t jY 6 / is survival as estimated by the uncorrected rescaling method, which inverts TO .s/. The corrected scaled estimator labelled ‘SP scaling’ in Figure 4 is the CIF GO 0 .t / D PO .Y 6 t / as given by equation (11). Our correction estimates survival only in the restricted domain t 6 . The restricted domain points to the inability of nonparametric methods to estimate the MGF or transmittance without additional assumptions [3, Theorem 4].

Appendix B. Transmittance calculations for linear hazard rate distribution extrapolation and proportional incidence With constant, independent of , and with extrapolation S.t / D S. / S0 .t / for t > utilising the LHRD survival S0 .u/ D exp ˛ u u2 =.2ˇ 2 / ; for u > 0, with ˛ > 0; ˇ > 0, equation (10), becomes Z T .s/ D S. / exp.s /

1

e

exp.s u/ dS0 .u/

0

D S. / exp.s / Œ1 C s

.ˇ.s ˛//

(12)

where ˇ > 0 and .u/ D

ˆ.u/ .u/


Statist. Med. 2014, 33 1621–1645

1643

Equation (12) uses the MGF of the LHRD distribution, the expression within braces on the line above, which can be calculated by completing the square of the exponent of the survival function S0 ; see [1] for

H. M. HUDSON ET AL.

its full derivation. Here, ˆ.u/ and .u/ are the distribution function and density of the standard normal distribution. Noting expressions for the derivatives 0 00

.u/ D 1 C u .u/; .u/ D .u/ C u 0 .u/

and setting r D ˇ.s ˛/, we have T e .s/ D S. / m.s/ with m.s/ D exp.s / Œ1 C ˇs

.r/ m .s/ D m.s/ exp.s / ˇ .r/ C ˇ 2 s 0 .r/ m00 .s/ D 2 m0 .s/ 2 m.s/ C 2ˇ 2 0 .r/ C ˇ 3 s 00 .r/ 0

Appendix C. Derivatives of the cumulant generating function in the randomisation-stroke-death model Because death is the sole endpoint, first-passage transmittance to endpoint is the MGF of survival time, and saddlepoint inversion requires derivatives of the CGGF, the logarithm of the MGF. We have K.s/ D log.MY .s// 0 0 D log T123 .s/ C T13 .s/ 0 0 K 0 .s/ D T123 .s/ C T13 .s/ =MY .s/ and 00 00 00 0 K .s/ D T123 .s/ C T13 .s/ =MY .s/ K .s/2 It suffices to determine component transmittances and their derivatives. Noting that T123 .s/ D T12 .s/ T23 .s/ we have 0 0 0 D T12 T23 C T12 T23 ; T123 00 00 0 0 00 T123 D T12 T23 C 2T12 T23 C T12 T23

Thus, we can compute K and its first two derivatives from the component transmittances.

Acknowledgements The authors thank the participants and the LIPID trial investigators. Hudson acknowledges reimbursement of research expenses for attendance and presentation at ISCB 2010, Montpellier, by the NHMRC CTC and by The George Institute. Heritier acknowledges the support of NHMRC Program Grant No. 571281. We are also grateful to Prof. Ronald Butler for early advice on saddlepoint methods and access to his textbook. The authors thank the referees for insightful comments and suggestions.

References

1644

1. Lo S, Heritier S, Hudson M. Saddlepoint approximation for semi-Markov processes with application to a cardiovascular randomized study. Computational Statistics and Data Analysis 2009; 53(3):683–698. (Special Issue in Clinical Trials). DOI: 10.1016/j.csda.2008.09.003. 2. Putter H, Fiocco1 M, Geskus RB. Tutorial in biostatistics: competing risks and multistate models. Statistics in Medicine 2007; 26:2389–2430. DOI: 10.1002/sim.2712. 3. Butler RW, Bronson DA. Bootstrap confidence envelopes for sojourn distributions in multistate semi-Markov models with right censoring. Biometrika 2012; 99(4):959–972. DOI: 10.1093/biomet/ass036. 4. Butler RW, Huzurbazar AV. Stochastic network models for survival analysis. Journal of the American Statistical Association 1997; 92:245–257. 5. Butler RW, Bronson DA. Bootstrapping survival times in stochastic systems by using saddlepoint approximations. Journal of the Royal Statistical Society, Series B 2002; 64:31–49. DOI: 10.1111/1467-9868.00323. 6. Butler RW. Saddlepoint Approximations with Applications. Cambridge University Press: Cambridge, 2007. DOI: 10.1017/CBO9780511619083. 7. Lagakos SW, Sommer CJ, Zelen M. Semi-Markov models for partially censored data. Biometrika 1978; 65:311–317. 8. Dinse G, Larson MG. A note on semi-Markov models for partially censored data. Biometrika 1986; 73:379–86.


Statist. Med. 2014, 33 1621–1645

H. M. HUDSON ET AL. 9. Nicolaie MA, van Houwelingen HC, Putter H. Vertical modeling: a pattern mixture approach for competing risks modeling. Statistics in Medicine 2010; 29:1190–1205. DOI: 10.1002/sim.3844. 10. The LIPID Study Group. Prevention of cardiovascular events and death with pravastatin in patients with coronary heart disease and a broad range of initial cholesterol levels. New England Journal of Medicine 1998; 339:1349–1357. 11. Huzurbazar AV, Williams BJ. Incorporating covariates in flowgraph models: applications to recurrent event data. Technometrics 2010; 52(2):198–208. DOI: 10.1198/TECH.2010.08044. 12. Lugannani R, Rice SO. Saddlepoint approximation for the distribution of the sum of independent random variables. Advances in Applied Probability 1984; 12:475–490. 13. Frydman H, Szarek M. Estimation of overall survival in an ‘illness-death’ model with application to the vertical transmission of HIV-1. Statistics in Medicine 2009; 29:2045–2054. DOI: 10.1002/sim.3949. 14. Stute W. The bias of Kaplan-Meier integrals. Scandinavian Journal of Statistics 1994; 21(4):475–484. 15. Stute W. The statistical analysis of Kaplan-Meier integrals. In Analysis of Censored Data, Vol. 27 of IMS Lecture Notes – Monograph Series, Koul HL, Deshpande JV (eds). Institute of Mathematical Statistics: Hayward, CA, 1995; 231–254. 16. Suzukawa A. Asymptotic properties of Aalen-Johansen integrals for competing risks data. Journal of the Japan Statistical Society 2002; 32(1):77–93. 17. Kodlin D. A new response time distribution. Biometrics 1967; 23:227–239. 18. Bain LJ. Analysis for the linear failure-rate life-testing distribution. Technometrics 1974; 16(4):551–559. 19. Kalbfleisch JD, Prentice RL. The Statistical Analysis of Failure Time Data, 2nd edn. Wiley: New York, 2002. 20. Putter H, Sasako M, Hartgrink HH, van de Velde CJH, van Houwelingen JC. Long-term survival with non-proportional hazards: results from the Dutch Gastric Cancer Trial. Statistics in Medicine 2005; 24:2807–2821. DOI: 10.1002/sim.2143. 21. Andersen PK, Esbjerg S, Sorensen TIA. Multi-state models for bleeding episodes and mortality in liver cirrhosis. Statistics in Medicine 2000; 19:587–599. DOI: 10.1002/(SICI)1097-0258(20000229)19:43.0.CO;2-0. 22. Satten GA, Datta S. Marginal estimation for multistage models: waiting time distributions and competing risk analyses. Statistics in Medicine 2002; 21:3–19. DOI: 10.1002/sim.967. 23. Cook RJ, Major P. Multistate analysis of skeletal events in patients with bone metastases. Clinical Cancer Research 2006; 12(20 Suppl):6264s–6269s. DOI: 10.1158/1078-0432.CCR-06-0654. 24. Ouhbi B, Limnios N. Nonparametric estimation for semi-Markov processes based on its hazard rate functions. Statistical Inference for Stochastic Processes 1999; 2:151–173. 25. Lagakos SW. Time-to-event analyses for long-term treatments–the APPROVe trial. New England Journal of Medicine 2006; 355:113–117. DOI: 10.1056/NEJMp068137. 26. Aalen O, Johansen S. An empirical transition matrix for nonhomogeneous Markov chains based on censored observations. Scandinavian Journal of Statistics 1978; 5:141–150. 27. Andersen PK, Borgan O, Gill RD, Keiding N. Statistical Models Based on Counting Processes. Springer Verlag: New York, 1993. 28. de Wreede LC, Fiocco M, Putter H. mstate: an R package for the analysis of competing risks and multi-state models. Journal of Statistical Software 2011; 38(7):1–30. DOI: 10.1016/j.cmpb.2010.01.001. 29. The LIPID Study Group. Long-term effectiveness and safety of pravastatin in 9014 patients with coronary heart disease and average cholesterol concentrations: the lipid trial follow-up. Lancet 2002; 359(9315):1379–1387. DOI: 10.1016/S0140-6736(02)08351-4. 30. Baigent C, Keech A, Kearney P, Blackwell L, Buck G, Pollicino C, Kirby A, Sourjina T, Peto R, Collins R, Simes R. The cholesterol treatment trialists’ (CTT) collaborators. Efficacy and safety of cholesterol-lowering treatment: prospective meta-analysis of data from 90,056 participants in 14 randomised trials of statins. Lancet 2005; 366:1267–1278. DOI: 10.1016/S0140-6736(05)67394-1. 31. Efron B. The two sample problem with censored data. In Proceedings of the 5th Berkeley Symposium, Vol. 4. University of California Press: Berkeley, 1967; 831–853. 32. Field C, Ronchetti E. Small Sample Asymptotics, vol. 13 of Lecture Notes – Monograph Series. Institute of Mathematical Statistics: Hayward, CA, 1990. 33. McCullagh P. Does the moment-generating function characterize a distribution? The American Statistician 1994; 48(3):208. 34. Meier P, Karrison T, Chappel R, Xie H. The price of Kaplan–Meier. Journal of the American Statistical Association 2004; 99(467):890–896. Review Article, DOI: 10.1198/016214504000001259. 35. Klein JP, Moeschberger ML. Survival Analysis: Techniques for Censored and Truncated Data. Springer: New York, 2003. 36. Royston P, Parmar MKB. The use of restricted mean survival time to estimate the treatment effect in randomized clinical trials when the proportional hazards assumption is in doubt. Statistics in Medicine 2011; 30:2409–2421. DOI: 10.1002/sim.4274. 37. Gill RD. Censoring and Stochastic Processes. Mathematish Centrum: Amsterdam, 1980.

1645


Statist. Med. 2014, 33 1621–1645

Semiparametric transformation models for semicompeting survival data.

Semiparametric methods for center effect measures based on the ratio of survival functions.

Combined estimating equation approaches for semiparametric transformation models with length-biased survival data.

Variable selection in strong hierarchical semiparametric models for longitudinal data.

Semiparametric methods to contrast gap time survival functions: Application to repeat kidney transplantation.

Semiparametric models for malaria rapid diagnosis test result.

Robust analysis of semiparametric renewal process models.

Effects of Photo and Genotype-Based Misidentification Error on Estimates of Survival, Detection and State Transition using Multistate Survival Models.

A comparison of methods to adjust for continuous covariates in the analysis of randomised trials.

Methods of designing two-stage winner trials with survival outcomes.

Comparison of semiparametric receiver operating characteristic models on observer data.

Causal estimation using semiparametric transformation models under prevalent sampling.

Semiparametric Bayes hierarchical models with mean and variance constraints.

Randomised trials for the Fitbit generation.

Semiparametric transformation models for causal inference in time to event studies with all-or-nothing compliance.

Randomised trials in developing countries.

Spatio-temporal semiparametric models for NO₂ and PM₁₀ concentration levels in Athens, Greece.

Multiple Testing under Dependence via Semiparametric Graphical Models.

Semiparametric analysis of linear transformation models with covariate measurement errors.

A Procedure for Deriving Formulas to Convert Transition Rates to Probabilities for Multistate Markov Models.

Randomised clinical trials.

Randomised consent trials.

Randomised controlled clinical trials.

Randomised controlled trials?