PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Stat Methods Med Res. Author manuscript; available in PMC 2017 March 24.
Published before final editing as:
PMCID: PMC5311029
NIHMSID: NIHMS758604

An ensemble survival model for estimating relative residual longevity following stroke: Application to mortality data in the chronic dialysis population

Abstract

Time-dependent covariates can be modeled within the Cox regression framework and can allow both proportional and nonproportional hazards for the risk factor of research interest. However, in many areas of health services research, interest centers on being able to estimate residual longevity after the occurrence of a particular event such as stroke. The survival trajectory of patients experiencing a stroke can be potentially influenced by stroke type (hemorrhagic or ischemic), time of the stroke (relative to time zero), time since the stroke occurred, or a combination of these factors. In such situations, researchers are more interested in estimating lifetime lost due to stroke rather than merely estimating the relative hazard due to stroke. To achieve this, we propose an ensemble approach using the generalized gamma distribution by means of a semi-Markov type model with an additive hazards extension. Our modeling framework allows stroke as a time-dependent covariate to affect all three parameters (location, scale, and shape) of the generalized gamma distribution. Using the concept of relative times, we answer the research question by estimating residual life lost due to ischemic and hemorrhagic stroke in the chronic dialysis population.

Keywords: Additive hazards, generalized gamma, relative times, residual median life, semi-Markov model

1 Introduction

Of the standard approaches available for the analysis of right censored survival data in biomedical applications, the semi-parametric Cox proportional hazards model1 is the most popular. This is partially due to the fact that in such applications research interest often centers on being able to estimate the relative risk of an event (say death) for patients belonging to one group as compared to patients belonging to another group. Here, group membership can be attributed to the type of treatment received (e.g. treatment versus placebo) or in many cases to the presence or absence of a particular demographic (gender, race, BMI group) or biological condition (cancer versus no cancer). The main strength of this approach is that through maximization of a partial likelihood, the estimated regression coefficients can be exponentiated to accord direct clinical interpretation of hazards relative to group membership irrespective of how hazards in each group change as a function of time. Proportionality of hazards also implies that the survival functions in the groups are related to each other by means of a power function, and this allows plotting survival profiles for patients belonging to each group when a baseline survival function is modeled or assumed. The validity of the proportional hazards assumption for covariates of interest can be assessed using standard approaches discussed in literature.2 When this assumption is violated, clinical interpretation of relative hazard is more involved than in the former case and alternative methods of analyses may be preferred.

Parametric methods of survival analyses permit the flexibility of modeling a variety of hazard shapes that can be described once the parameters specifying the particular distribution have been estimated by maximization of the overall likelihood function. With the advancement of statistical software, many parametric distributions can be used to model survival data thereby facilitating between-group comparison for individuals exposed to a condition as opposed to individuals not exposed to that condition through the concept of relative times (RT). Specifically, RT is the ratio of times that a given proportion (or percent) of individuals belonging to one group experience an event of interest as compared to individuals not belonging to that group. Parametric models, therefore, can have a bigger impact in those biomedical applications where comparative estimation of event-free survival times for individuals exposed to a medical condition (say group-1) versus individuals not exposed to this condition (group-2) is the focus of research interest. For example, in certain cancer studies researchers may be more interested in estimating the residual longevity of patients possessing certain characteristics compared to patients who do not possess those characteristics. The hazards in patients with and without these characteristics may change differentially across time making the quantification of such comparative estimation difficult—even with the allowance of time-dependent covariates in a Cox regression framework. In such applications, parametric models can be deployed to answer the research question provided reasonably accurate inference can be drawn about the various hazard shapes attributable to the estimated parameters of these distributions. Additionally, hazard functions for the two groups of individuals can also be plotted over time by deriving these functions from the parametric survival function.

Our work is primarily motivated by research interest in estimating residual longevity in the chronic dialysis population due to two different types of stroke: hemorrhagic and ischemic. As the two types of strokes are physiologically different, it is hypothesized that they will confer completely different survival experience—both from each other and from those not experiencing either event. Here, stroke is a time-dependent covariate, that is, dialysis patients who enter the cohort do so with a “no stroke” status that may or may not change to one of the two “stroke” (hemorrhagic or ischemic) statuses at some point in the observation window. Note that some subjects may have incurred a stroke prior to dialysis initiation, but as Medicaid data is not available on all patients prior to dialysis initiation, we have focused our interest in studying the impact of ischemic and hemorrhagic strokes that occur while subjects are on dialysis. A priori, it is not known whether stroke affects patient survival through its time-dependent occurrence alone or its effect is accentuated by time spent since the stroke occurred. From a clinical perspective, doctors would like to inform patients as to the expected median residual life lost (RLL) due to the occurrence of each type of stroke, allowing also for adjustment for a host of demographic risk factors and comorbidities. Simultaneously, researchers would prefer to retain conventional hazard ratio (HR) interpretations for these accompanying risk factors. Another research interest is to study and compare the hazard ratios, which are potentially nonconstant of the two types of strokes over time so as to evaluate the differential relative risks that the two stroke types accord over time.

In order to answer the research question, we adopt an ensemble approach that combines various aspects of standard and advanced survival models. Section 2 summarizes information related to the motivating example. Section 3 provides model building justification for our model. In section 4, we briefly discuss the three-parameter generalized gamma (GG) distribution and its properties in the context of our example. Main results related to RT calculations are discussed in section 5. In section 6, we summarize the results and discuss potential extensions to our approach.

2 Motivating example

In this section, we briefly summarize background information related to this project; details of cohort construction are given in Wetmore et al.3 Briefly, a large cohort with 69,371 patients representing the chronic dialysis population in the time period 2000–2005 who were eligible for receiving both Medicare and Medicaid was constructed by linking the United States Renal Data System4 (USRDS) to Medicaid claims data. The research team developed clinical algorithms that allowed identification of incident hemorrhagic and ischemic strokes from Medicare claims data.5,6 Follow-up for this cohort began upon dialysis initiation (plus 90 days, as is standard with USRDS data). Subjects were followed until death, the outcome of interest, and subjects were right censored in the case they lost their Medicare or Medicaid coverage, received a kidney transplant, or were otherwise lost to follow-up (follow-up ended on 31 December 2005). Stroke and atrial fibrillation were treated as time-dependent covariates whereas all other risk factors were treated as time-independent baseline covariates. Identification of stroke type (hemorrhagic or ischemic) was done using international classification of diseases-9 (ICD-9) codes, and for those with more than one stroke event, only the first event was used. This approach allowed identification of 2391 ischemic and 534 hemorrhagic strokes. As the main research question concerns estimation of residual lifetime lost due to the two stroke types, all subsequent section results will be discussed mainly with respect to stroke as a time-dependent covariate; the adjustment for other risk factors only being discussed in the model building part of sections 3 and 4.

3 Semi-Markov model with additive hazard extension

3.1 Initial model building assessments

This section focusses on the model building procedure with reference to using a semi-Markov model with an additive hazard extension for our problem. As mentioned in section 1, a priori it is not known whether relative risk of death due to the occurrence of a stroke is affected by time at which stroke occurs. A Kaplan–Meier (KM) curve (see Figure 1) shows the survival profile of patients suffering hemorrhagic and ischemic stroke after occurrence of the corresponding stroke type (i.e. with time axis representing ”time since stroke”). When stratified by year on dialysis of stroke occurrence (see Figure 2(a) for hemorrhagic, Figure 2(b) for ischemic), the corresponding log-rank tests suggest that stroke vintage is not associated with time to death following a stroke (p-values are 0.530 and 0.401, respectively). This further suggests that the predicted survival profiles following a stroke—regardless of when the stroke occurred—from our final model should somewhat resemble the shape of the KM curves in Figure 1 when the values of all covariates are adjusted to match the overall profile of the patients that experienced strokes in the cohort. The goodness of fit for any intermediate model choice can, therefore, be assessed diagnostically in terms of visual comparison of model-based predicted survival curves to the corresponding KM curve in Figure 1. Such an approach of diagnostic assessment in the time-independent covariate scenario has been adopted by Storer et al.7 for their Cox PH type model. To summarize to this point, we assumed time of stroke occurrence is not associated with subsequent mortality, and have identified “observed” plots to utilize for assessment of subsequent models built.

Figure 1
Kaplan–Meier curve illustrating survival after an ischemic or hemorrhagic stroke.
Figure 2
Kaplan–Meier curve stratified by year of occurrence for survival after first stroke: (a) hemorrhagic, (b) ischemic.

As a next step, we deployed the standard Cox PH model with stroke and atrial fibrillation incorporated as time-dependent covariates. Without loss of generality we describe how stroke was built into this phase of the modeling by including only parameters for stroke using

h(tZ)=h0(t)exp{β1Zt+β2Ztf(t)+β3X1+β4X2++βkXk2}
(1)

Here, X1, X2, …, Xk−2 are the other risk factors including atrial fibrillation, β3, β4, …, βk are the regression coefficients associated with these other risk factors, Zt is the time-dependent indicator variable for stroke (1 if stroke occurs, else 0), β1 is the regression coefficient associated with Zt, whereas Zt · f(t) is the time-dependent term representing a functional interaction form between stroke and time, and β2 is the regression coefficient associated with Zt · f(t). Note that we have used the subscript t on Z to indicate the time-dependency of this variable. As we have assumed that time at which stroke occurs does not affect survival, we can replace Zt · f(t) by Zt · f(t − s) where s = time at which stroke occurs implying that post stroke, survival is affected only by: (1) the fact that stroke has occurred (Zt) and (2) the time spent since occurrence of stroke Zt · f(t − s).

For an indicator variable Zt (equals 1 once stroke occurs, 0 otherwise), different functional forms of Zt · f(t − s) can be modeled with Zt · f(t − s) = Zt · (t − s) and Zt · log(ts) being the most common. Hougaard8 has discussed how in the former case for s = 0, h(t|Z) = h0(t)exp(β1Zt + β2Zt · t) yields a relative risk function that has the same form as that of a Gompertz parametric model with different values of the shape parameter for the two groups (stroke versus no stroke). Likewise, in the latter case for s = 0, h(t|Z) = h0(t)exp{β1Zt + β2Zt · log(t)} yields a relative risk function that has the same form as that of a Weibull model with different values of the shape parameter for the two groups (stroke versus no stroke). Substituting s = 0 and generating a predicted survival profile for patients who experience a stroke (this is done for both stroke types), we find that both functional forms prove inadequate in passing the graphical diagnostic check mentioned above. This implies that although stroke as a time-dependent covariate can be incorporated in the Cox regression framework, it is difficult to find the appropriate functional form that satisfies the diagnostic checking procedure. Thus, testing for the proportionality of hazard using the Therneau—Grambsch nonproportionality test9 alone is not enough as a significant result for this test might be indicative of a number of possible model failures. As noted by Keele,10 the correct interpretation of a significant test result on the nonproportionality test is not evidence that the hazards are not proportional but instead that if the covariate-dependent specification is correct, then there is evidence that hazards are not proportional. Thus, when the model is misspecified, the Therneau–Grambsch test may yield a false positive result.11,12 The latter paper also discusses importance of incorporating possible nonlinear covariate effects as well as considering other parametric survival models that result in nonproportional hazards. In our paper (in the next section), we therefore, use the GG family of distributions to model the nonproportional hazards part of our final model.

3.2 A state transition model building approach

One way to handle this problem is to move from a time-dependent covariate approach to an equivalent state-dependent transition approach. That is, we exploit the fact that a single time-dependent covariate (in this case stroke) can also be expressed as a single intermediate state transition model as represented in Figure 3. With reference to this figure, State-1 is the initial state representing entry into cohort, State-2 is the single intermediate transition state representing stroke occurrence, and State-3 is the single absorbing state signaling death. Entry into cohort takes place at time t = 0 from which an individual can progress to death at time t either directly (those with no strokes) or through the single transitional State-2 (a stroke) occurring at time s such that ts. Individuals who do not experience the death are considered right censored at time t. Under the Markov assumption which states that what happens to an individual after he has visited a particular state depends on their current state and not on the preceding history, it can be shown13 that the transition hazard hij(t) for moving from state i to state j can be modeled as

hi3(t)=h30(t)exp{β13X+(i1)ΔX+(i1)δ}i=1,2
(2)

Here, h30(t) represents some baseline hazard attributable to moving into State-3, β13 represents the vector of regression coefficients for the other risk factors, X is the vector of other risk factors, Δ represents the difference in covariate effects of these other risk factors for death after having stroke, compared to before having stroke, and δ represents the estimate of the stroke effect under the Markov assumption. This stringent Markov assumption, however, can be relaxed to form a semi-Markov model by allowing the sojourn times (incorporated through covariates) to depend on the history of the process through: (i) the current state itself and (ii) time spent in that state. Thus a semi-Markov model allows us the flexibility to incorporate the influence of time spent since stroke on the hazard of death for those that have had a stroke. It should be noted that though a semi-Markov model can be further extended to form an extended semi-Markov model (where in addition to the semi-Markov assumptions, entry time in a particular state can also affect hazard of death), such an extension is unnecessary in our case due to the assumption that the time at which stroke occurs does not play a role toward affecting patient survival from the time of the stroke until death.

Figure 3
Single-state transition model for mortality with stroke as an intermediate state.

In the context of Figure 3 and the fact that stroke is a time-dependent (and equivalently state-dependent) covariate, meaning the contribution of (ts) to the hazard of death is valid only for those had strokes (as those that have not had them do not pass through State-2), the semi-Markov model can be augmented with an additive hazards component in the following way

hi3(t)=h130{min(t,s)}exp(β13X)+h23(ts)I(ts)i=1,2
(3)

where s is the instant just before the occurrence of a stroke at time = s and I is the indicator function that equals 1 if ts and 0 otherwise.

Thus for those with no strokes, i = 1 and t < s resulting in

h13(t)=h130(t)exp(β13X)
(4)

and for those with strokes, i = 2 and ts resulting in

h23(t)=h130(s)exp(β13X)+h23(ts)
(5)

Here, h130(t) is some baseline hazard associated with the transition from State-1 directly to State-3. However, for those who experience a stroke at time s, the corresponding part of the baseline hazard is only h130(s) because after the stroke they experience a change in their overall hazard owing to the effect of stroke, measured as a function of time spent since stroke. This incremental hazard is given by h23(ts) and it can be modeled in two ways (depending on how the covariates affect this hazard function), the discussion of which is postponed to section 4.4. In the next section, we show how we can use the three parameter GG distribution to identify h130(s) and h23(t − s) and how the properties of this distribution can be used to do the calculations related to RLL due to each stroke type using the concept of RT. More examples of multistate models with additive hazards extensions can be found in Klein et al.14 and in Shu et al.15 Variance estimation for some of these models can be found in Klein et al.16 Large-sample properties of three different additive hazards extension models can be found in Shu.17 Likewise, examples with a nonparametric semi-Markov model and a semi-Markov model incorporating covariate information are illustrated in Voelkel et al.18 and Andersen et al.,19 respectively. Estimation and prediction in a semi-Markov model in a Cox regression framework with baseline hazards and covariates depending on recurrence times has been illustrated in Dabrowska.20

4 Using the GG distribution for calculating residual life lost due to stroke

4.1 The generalized gamma distribution as a family of distributions

The GG distribution21 is a three-parameter distribution family with location parameter μ, scale parameter σ, and shape parameter λ that generalizes the two-parameter gamma distribution. The general notation for specifying the distribution is GG(μ, σ, λ). Its density function is

fGG(t)=λσtΓ(λ2)[λ2{exp(μ)t}λσ]λ2exp[λ2{exp(μ)t}λσ]
(6)

where σ > 0, μ [set membership] (−∞, ∞), λ [set membership] (−∞, ∞), and Γ(x)=0mx1em dm is the gamma function of x.

A complete taxonomy of the various hazard functions for the GG family has been explained in Cox et al.22 and only the relevant aspects are discussed here. Briefly, the GG family allows the flexibility of modeling different shapes of hazard functions such as increasing from 0 to ∞, increasing from a constant to ∞, decreasing from ∞ to 0, decreasing from ∞ to a constant, arc-shaped hazards, and bathtub-shaped hazards. Also many popular parametric distributions are special case members of this family. Thus, λ = σ gives a two-parameter gamma distribution; μ = 0; σ = 1 gives a standard gamma distribution for fixed λ; λ = 1 gives a Weibull distribution; λ = σ = 1 gives an exponential distribution; the limiting case λ = 0 gives a lognormal distribution; λ = −1 gives an inverse Weibull distribution; λ = −σ gives an inverse gamma distribution; λ = 1/σ gives an ammag distribution; λ = −1/σ gives an inverse ammag distribution; and a lognormal distribution with σ′ = 1:82σ approximates the loglogistic distribution.22 Maximum likelihood estimation using standard software can be used to obtain estimates of the GG distribution parameters with parsimonious reductions resulting in fitting of well-known distributions.

4.2 Quantiles of the GG distribution

As the main research interest is to quantify residual longevity, it is informative to study the quantiles of a GG(μ, σ, λ) distribution. Cox et al.22 show that

log{tGG(μ,σ,λ)(p)}=μ+σlog{tGG(0,1,λ)(p)}=μ+σgλ(p)
(7)

Here, gλ(p) is the logarithm of the 100pth quantile (0 ≤ p ≤ 1) from the GG(0, 1, λ) distribution. The location parameter μ acts as a time-multiplier and governs the values of the median for fixed values of σ and λ, resulting in an accelerated failure time (AFT) model. The scale parameter σ determines the interquartile ratio for fixed values of λ and independently of μ. The shape parameter λ determines the GG(0, 1, λ) distribution. Together, σ and λ describe the type of hazard function for the GG(μ, σ, λ) distribution. Operational details of a standard AFT model can be found in Wei.23

4.3 Relative times

The time by which 100 p% of the population experience an event can lead to a statistic called “relative times RT(p)”, which can be used to compare survival profiles of patients in different groups (e.g. those that had a stroke versus those that did not). Thus

RT(p)=t1(p)t0(p)=S11(1p)S01(1p)
(8)

where Si1(1p) is the inverse survival function for group i (i = 0,1).

The interpretation of RT(p) is that the time required for 100p% of individuals in the stroke group to experience death is RT(p) times the time required for 100p% of individuals in the no stroke group to experience death. Thus if (μ0, σ0, λ0) and (μ1, σ1, λ1) denote two different sets of GG parameter values, then

RT(p)=exp{(μ1μ0)+σ1gλ1(p)σ0gλ0(p)}
(9)

The manner in which covariates affect RT(p) can be summarized as:

  • (1)
    If λ1 = λ0 and σ1 = σ0 then we have a conventional AFT model resulting in nonproportional hazards, but proportional RT; that is, covariates affect μ only.
  • (2)
    If only λ1 = λ0 then we have a model that results in nonproportional hazards and nonproportional RT(p); i.e. covariates affect both μ and σ.
  • (3)
    Full generalization is obtained by having covariates affect all three parameters.
  • (4)
    Reduced parsimonious models result in the fitting of family members of the GG distribution.

4.4 Step-wise approach to identify parameters for the motivating example

With reference to the model in (3), following step-wise approach can be utilized to answer the research question.

  • (1)
    For pre-stroke times fit the Cox PH model given in (4) with right-censoring also occurring at the time of stroke (when applicable). Obtain estimates of β13 for all risk factors other than stroke. The Cox PH framework also allows estimation of risk due to AF, which is a time-dependent covariate.
  • (2)
    Using maximum likelihood estimation method, parametrically estimate the parameters for the most parsimonious subfamily of GG distribution for these pre-stroke times (including subjects with no strokes prior to death) in the absence of any covariates.
  • (3)
    Then generate baseline survival curve for pre-stroke times using the general relation
    S0(t)=exp{0th0(t)dt}using either step 1 or step 2.
  • (4)
    Generate a survival profile for pre-stroke times using the average values of all other covariates X, using the (general) relation S0,x(t)={S0(t)}β13X
  • (5)
    For any time s (of interest) for occurrence of a stroke, calculate h130(s). Adjust this for the effect of all other covariates X.
  • (6)
    The hazard h23(ts) for those that had strokes as a function of “time since stroke” can be obtained in two ways.
    • (a)
      The first way is to have all the other covariates affect μ in the AFT model of (7) and then estimate parameters of the most parsimonious subfamily of the GG distribution. Software (e.g. the SAS procedure LIFEREG) can be used for this purpose. Then generate survival curves for those that had strokes as a function of ts, and augment them to the survival curve for those that had not had a stroke at time of stroke = s. These survival curves can be generated for either particular values of covariates (such as African American males ≥ 65 years, say) or for an average profile of patients in the cohort as done in step 4 of this list.
    • (b)
      The second way is to consider h23(ts) = h230(ts) · exp(β13X + Δ′X) such that only the baseline hazard h230(ts) as a function of “time since stroke” is due to the GG distribution and Δ is the increment vector of regression coefficients β13 corresponding to the other covariates that indicates whether or not the effect of these risk factors on mortality change due to the occurrence of a stroke. As β13X now occurs on both sides of the + sign in (5), software such as the SAS procedure NLMIXED can be used to maximize the overall likelihood (see Cox et al.22 for sample code using NLMIXED) resulting from (5). Then, similar to step 6a above, the desired survival profiles can be generated. However, it should be noted that the PH assumption needs to be validated here.
  • (7)
    Residual longevity following a stroke can be readily calculated in terms of the 100pth percentile (for median, p = 0.5) survival time after the occurrence of a stroke at time s from step 6 (a or b) above. With reference to Figure 4, this corresponds to calculating
    tstroke(p)s=Sstroke1{(1p)Snostroke(s)}s.
    Figure 4
    Calculating residual longevity following stroke occurrence at time = s.
  • (8)
    Residual longevity absent any strokes for any time s can be calculated by evaluating the 100pth percentile of survival time from the conditional survival function S0(t|ts). With reference to Figure 4, this corresponds to calculating
    tnostroke(p)s=Snostroke1{(1p)Snostroke(s)}s.
  • (9)
    Then RT(p)=tnostroke(p)ststroke(p)s can be used to calculate the relative 100pth percentile of time required for absent a stroke to experience death as compared to when a stroke occurs. Alternatively, the 100pth percentile of RLL due to stroke, a concept more appealing in practice to the clinicians can be calculated using the relation given by RLL(p)=Snostroke1{(1p)Snostroke(s)}Sstroke1{(1p)Snostroke(s)}. Specifically, RLL(0.25), RLL(0.5), and RLL(0.75) may prove to be quite informative in explaining the skew in the RLL due to stroke. These calculations are graphed in Figure 4. It should also be noted that though Figure 4 depicts the comparison of median (assuming p = 0.5) residual lifetime of those having a stroke (at say time s) to the median residual lifetime of those not having a stroke under the assumption that they will always maintain their “no-stroke” status, it can also be used for alternative calculations. Thus acknowledging that some patients in the “no-stroke” state will later go on to have a stroke, we can also use it to compare the RLL between those that have a stroke (at time s) and those who are stroke-free at time s but may experience a stroke a future time s′. For example, say p=0.5, s=12, and s′ = 36, then the median RLL can now be calculated as
    Sstrokes=361{(0.5)Snostroke(12)}Sstrokes=121{(0.5)Snostroke(12)}.
  • (10)
    Evaluation of the model fit can be performed using the graphical diagnostic checking procedure mentioned in section 3.1.

In the above 10-step procedure, we have explained how RLL(p) and RT(p) can be estimated with all covariate values adjusted at their average values (representing average survival profiles of patients in the stroke and no-stroke groups). We can also generate survival profiles for any particular combination of covariate values (say, white male, age 65 with hypertension and diabetes) and therefore estimate RLL(p) and RT(p) for those combination of covariate values.

5 Results

Descriptive statistics for risk factors related to patient demographics can be found in Wetmore et al.24 Here, we discuss the HR estimates of these risk factors only in relevance to our model of (3). Table 1[A] displays the results representing (4) for no or pre-stroke times using a Cox PH model with AF treated as a time-dependent covariate. Baseline hazard generated from this model suggests a Weibull distribution as a parsimonious reduction of a GG distribution (λ = 1) for the times prior to stroke events, with a gradually decreasing hazard that slowly stabilizes over time with parameter estimates given by μ^ = 4.385, σ^ = 1.096 (obtained from a linear fit of logarithm of the cumulative hazard versus logarithm of time). The survival curve shown in Figure 4 was generated using these parameters by adjusting the baseline survival curve for an average risk profile of patients in the cohort.

Table 1
Modeling results using equation (4) for [A] no or pre-stroke times and equation (5) for times following [B] hemorrhagic and [C] ischemic strokes.

Table 1[B] displays the results corresponding to the additive part of the hazard outlined in (5) used to model the survival experience following a hemorrhagic stroke (N = 534 such events) using the first approach (step 6a in section 4.4), and other risk factors affect only μ in the AFT model of (7) using PROC LIFEREG). From this table, it can be seen that a GG distribution with parameters estimates given by μ^ = −2.179 (99% CI: [−4.761, 0.403]), σ^ =1.258 (99% CI: [1.065, 1.485]), and λ^ = −2.281 (99% CI: [−2.853, −1.710]) provides the best fit (AIC 1681.89). Adjusting μ^ to the average risk profile of patients allows generation of Figure 4 for hemorrhagic strokes (this adjusts the location estimate to μ^ = −1.205). We also fitted the PH model mentioned in step 6(b) in section 4.4 using PROC NLMIXED. On doing so, we found that only two risk factors belonging to the Δ vector (age and smoking status were statistically significant at the 0.01 level). However, the Cox PH assumption did not hold true for these two variables and hence this model was discarded. (Note that this model requires estimation of 44 regression coefficients in addition to three parameters each for the nonstroke and hemorrhagic stroke groups which may be unreliable due to there being only N=409 death events in hemorrhagic stroke group).

Table 1[C] displays the results corresponding to the additive part of the hazard outlined in (5) used to model the survival experience following ischemic strokes (N=2381) using the first approach (other risk factors affect only μ in the AFT model of (7) using PROC LIFEREG). A lognormal distribution with parameters estimates given by μ^ = 3.041 (99% CI: [1.665, 4.416]) and σ^ =1.883 (99% CI: [1.793, 1.979]) provides the best fit (AIC 7411.36 for lognormal as compared to AIC 7412.29 for GG). Adjusting μ^ to the average risk profile of patients allows Figure 4 ischemic stroke results (this adjusts the location estimate to μ^ = 2.333). Here again, we fitted the PH model mentioned in response numbered 6b in section 4.4 using PROC NLMIXED and found that only two risk factors belonging to the Δ vector (ambulation and employment status were statistically significant at the 0.01 level). However, as the Cox PH assumption did not hold true for these two variables as evidenced by log—log survival plots, this model was discarded. It should be noted that though in our application the model using NLMIXED was not very useful, it can turn out to be the model of choice in other applications.

Tables 2 and and33 show calculations related to RT(p) and RLL(p) for different values of p when the two stroke types are separately compared to no stroke. The 99% CIs are generated using the delta method. All comparisons are adjusted for the average demographic and risk profile of patients in the cohort. From this output, we see that for patients who experienced a hemorrhagic stroke right at the entry point into the cohort, RLL(p) ranges from 2.291 months (5th percentile) to 53.453 months (60th percentile) with an estimated value of 39.903 months for median residual life lost. Likewise, for patients who experience a hemorrhagic stroke at 12 months into the cohort, the RLL(p) ranges from 2.919 months (5th percentile) to 42.461 months for median residual life lost. As approximately 38% of the no stroke group do not die in the 60 month observation window, RLL(p) calculations are restricted to the 60th percentile for s=0, 50th percentile for s=12, and 40th percentile for s=24. The increase in RLL(p) values for different values of s=0, 12, 24 is due the fact that though the overall survival experience following a stroke remains the same, the Weibull hazard of the reference (no stroke) group decreases gradually over time. Thus there is greater loss in residual life for patients who experience a stroke late into the observation window (time on dialysis) as opposed to early. Analogously, for patients who experienced an ischemic stroke right at the entry point into the cohort, RLL(p) ranges from 1.946 months (5th percentile) to 38.975 months (60th percentile) with an estimated value of 30.715 months for median residual life lost. Likewise, for patients who experience an ischemic stroke at 12 months into the cohort, the RLL(p) ranges from 2.575 months (5th percentile) to 33.272 months for median residual life lost. From these numbers (and from Figure 4), we can conclude that a hemorrhagic stroke confers a higher disadvantage in terms of residual life lost than an ischemic stroke. For example, for s=12, median RLL due to a hemorrhagic stroke exceeds median RLL due to an ischemic stroke by approximately nine months.

Table 2
Comparison of hemorrhagic stroke versus no stroke using RT(p) and RLLmonths(p).
Table 3
Comparison of ischemic stroke versus no stroke using RT(p) and RLLmonths(p).

The RLL(p) estimates are complemented by the RT(p) estimates in that they provide an interpretation of the multiplicative factor that the residual longevity would have increased for the 100pth percentile of patients who experience a stroke at time s had they continued on dialysis absent that stroke. Thus, for patients who experience a hemorrhagic stroke at s=12, RT(p) increases from 25.204 for p=0.05 to 54.025 for p=0.25 before gradually dropping down to 39.271 at p=0.5. Analogously for s=12, RT(p) for those that had ischemic stroke increases from 6.538 at p=0.05 to 6.814 at p=0.1 before decreasing gradually to 4.231 at p=0.5 suggesting that hemorrhagic stroke is a more drastic event than ischemic stroke. These estimates increase for patients who experience a stroke later in the observation window. Keeping p fixed, the confidence intervals are wider for both RT(p) and RLL(p) for high values of s (smaller risk set) as compared to small values of s (larger risk set).

6 Discussion

By using a semi-Markov model with an additive hazard extension, we were able to estimate relative residual longevity due to the two types of strokes. Estimates of RT(p) and RLL(p) obtained from this framework suggest that both stroke types are major, watershed events that radically change the survival profile of those who experience them. In fact, these stroke events affected all three parameters of the GG distribution model, which was then used to estimate the survival trajectory of patients who experienced strokes, which is essential to calculating residual longevity. This would not have been possible in a standard Cox PH model with stroke treated as a time-dependent covariate because in that case stroke would affect only the magnitude of the hazard in a proportional manner without affecting the shape. Even if the general Cox model framework can incorporate an interaction term of stroke with some functional form of time to account for the nonproportionality of hazard, it would be very difficult to reliably find this functional form. On the other hand, proportionality of hazard functions in the GG family is only limited to the special case of conventional Weibull regression models where the covariate of interest affects only the location parameter. Instead when the covariate can affect two or all three parameters of the GG family, our modeling framework allows the flexibility of using commonly available standard statistical software to model the nonproportionality of hazards and/or the nonproportionality of relative times in order to generate the survival profiles of interest (for any combination of covariate values).

Though we have focused on the nonproportional RT(p) and RLL(p) in this paper, it is also possible to estimate the nonproportional hazards due to stroke, and this has been done in Wetmore et al.24 as clinicians are more accustomed to the hazard ratio interpretation of risk factors. It is the motivating example, however, that has driven our current approach in the sense that it makes it possible for researchers to inform doctors and patients who experienced a stroke while on dialysis, the estimated median (or any other 100pth percentile) residual life lost due to the stroke. This is a paradigm shift from the usual approach in which only the mortality hazard ratio is reported.

Another interesting feature of our model is revealed by observing the exponentiated location estimates in Table 1[B] and [C]. The drastic nature of a hemorrhagic stroke is revealed by the fact that none of the other risk factors play a significant role in altering the risk of mortality in its presence. This is understandable given that median time to death after hemorrhagic stroke is only 0.84 months. On the other hand, in case of ischemic stroke, advancing age and white race still continue to significantly alter the mortality risk. In the context of the AFT extension of (5), every one year increase in age decreases the median time to death by a factor of 0.96, and race (black versus white) increases the median time to death by a factor of 1.75. A similar effect of race and age is observed in the no or pre-stroke time, and this can be verified using the relation between parameter estimates from a Cox PH model and an AFT Weibull model for no/pre-stroke survival given by β^PH = β^AFT−Weibull/σ. This suggests that after an ischemic stroke, race and age continue to retain their effect on mortality whereas the other risk factors are shadowed by the drastic nature of an ischemic stroke.

One limitation of our study is that the observation window is limited to a maximum possible follow-up time of 69 months resulting in a Weibull distribution with a gradually decreasing hazard for the nonstroke subjects. A larger follow-up window would result in the hazard for nonstroke subjects increasing gradually over time owing to the fact that this cohort represents a chronic dialysis population. In that case, a Weibull distribution perhaps, would not be the best choice and the new choice would affect the RT(p) and RLL(p) calculations accordingly. It is also possible that in that case we would have to consider the complicated extended semi-Markov model with time at which stoke occurs also playing a key role in the subsequent calculations. Another minor limitation of our study is that it uses claims-based data (all dates have been calculated using clinically approved algorithms rather than the exact event dates).

Our modeling framework also provides the flexibility to extend our approach to analysis for chronic dialysis cohort observed after 2005 and treated with exposure to new treatments (say, warfarin) with current cohort survivors left truncated at time of entry into the “new treatment-exposure” cohort. In this context, it may be interesting to investigate the effect of new treatment by studying RT(p), RLL(p), and nonproportional hazards across different time periods. Perhaps, the more involved model using PROC NLMIXED (discussed in section 4.4) could provide more useful in this regard. We envisage that our ensemble modeling framework will provide real life interpretations of research questions where the primary focus lies in drawing meaningful conclusions related to the survival and hazard trajectory of cohort participants using multi-parameter distributions. To this end, our model provides a reasonable solution.

Acknowledgments

Funding

This work was supported by National Institute of Diabetes and Digestive and Kidney Diseases (R01 DK080111); the National Kidney Foundation Young Investigator Award and the Sandra A. Daugherty Foundation Grant (K23 DK085378-01).

Footnotes

Conflict of interest

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

1. Cox DR. Regression models and life tables (with discussion) J Roy Stat Soc B. 1972;34:187–220.
2. Klein JP, Moeschberger ML. Survival analysis - Techniques for censored and truncated data. 2nd ed. Springer; New York: 2003. p. 304.
3. Wetmore JB, Ellerbeck EF, Mahnken JD, et al. Atrial fibrillation and risk of stroke in dialysis patients. Ann Epidemiol. 2013;23:112–118. [PMC free article] [PubMed]
4. United States Renal Data System . USRDS 2006 Annual Data Report: Atlas of end-stage renal disease in the United States. National Institutes of Health, National Institute of Diabetes and Digestive and Kidney Disease; 2006.
5. Wetmore JB, Ellerbeck EF, Mahnken JD, et al. Stroke and the ‘stroke belt’ in dialysis: contribution of patient characteristics to ischemic stroke rate and its geographic variation. J Am Soc Nephrol. 2013;24:2053–2061. [PubMed]
6. Wetmore JB, Phadnis MA, Mahnken JD, et al. Race, ethnicity, and state-by-state geographic variation in hemorrhagic stroke in dialysis patients. Clin J Am Soc Nephrol. 2014;9:756–763. [PubMed]
7. Storer BE, Gooley TA, Jones MP. Adjusted estimates for time-to-event endpoints. Lifetime Data Anal. 2008;14:484–495. [PMC free article] [PubMed]
8. Hougaard P. Analysis of multivariate survival data. Springer; New York: 2001. p. 87.
9. Grambsch PM, Therneau TM. Proportional hazards tests and diagnostics based on weighted residuals. Biometrika. 1994;81:515–526.
10. Keele L. Proportionality difficult: Testing for nonproportional hazards in Cox models. Political Anal. 2010;18:189–205.
11. Therneau TM, Grambsch PM, Fleming TR. Martingale based residuals for survival models. Biometrika. 1990;77:147–160.
12. Therneau TM, Grambsch PM. Modeling survival data: Extending the Cox model. Springer Verlag; New York: 2000.
13. Putter H, Fiocco M, Geskus RB. Tutorial in biostatistics: Competing risks and multi-state models. Stat Med. 2007;26:2389–2430. [PubMed]
14. Klein JP, Shu YY. Multi-state models for bone marrow transplantation studies. Stat Meth Med Res. 2002;11:117–139.
15. Shu YY, Klein JP. Additive hazards Markov regression models illustrated with bone marrow transplant data. Biometrika. 2005;92:283–301.
16. Klein JP, Qian C. Modeling multistate survival illustrated in bone marrow transplantation; Proceedings of the Biometrics Section of the American Statistical Society; 1996.pp. 93–102.
17. Shu YY. Multistate survival models: Theory and applications. Division of Biostatistics, The Medical College of Wisconsin; USA: 2001. PhD Dissertation.
18. Voelkel JG, Crowley J. Nonparametric inference for a class of semi-Markov processes with censored observations. Ann Stat. 1984;12:142–160.
19. Andersen PKA, Esbjerg S, Sorensen TIA. Multistate models for bleeding episodes and mortality in liver cirrhosis. Stat Med. 2000;19:587–599. [PubMed]
20. Dabrowska DM, Sun G, Horowitz MM. Cox regression in a Markov renewal model: An application of bone morrow transplant data. J Am Stat Assoc. 1994;89:867–877.
21. Stacy EW, Mihram GA. Parameter estimation for a generalized gamma distribution. Technometrics. 1965;7:349–358.
22. Cox C, Chu H, Schneider MF, et al. Parametric survival analysis and taxonomy of hazard functions for the generalized gamma distribution. Stat Med. 2007;26:4352–4374. [PubMed]
23. Wei LJ. The accelerated failure time model: A useful alternative to the Cox regression model in survival analysis. Stat Med. 1992;11:1871–1879. [PubMed]
24. Wetmore JB, Phadnis MA, Ellerbeck EF, et al. Stroke as a watershed event: Impact of first acute stroke on mortality in dialysis patients. Clin J Am Soc Nephrol. 2014;10:80–89. [PubMed]