Analysis of Heart Transplant Survival Data Using Generalized Additive Models

Abstract

The Stanford Heart Transplant data were collected to model survival in patients using penalized smoothing splines for covariates whose values change over the course of the study. The basic idea of the present study is to use a logistic regression model and a generalized additive model with 𝐵-splines to estimate the survival function. We model survival time as a function of patient covariates and transplant status and compare the results obtained using smoothing spline, partial logistic, Cox's proportional hazards, and piecewise exponential models.

1. Introduction

Cox’s proportional hazards model has been proposed based on the relationship between survival and the patient characteristics observed when the patient entered the study [1]. When the values of covariates change over the course of the study, however, a number of theoretical problems with respect to the baseline survival function and the baseline cumulative hazard function need to be solved [2]. Several prognostic models [3–6] have become as widely used as Cox’s proportional hazards model for the analysis of survival data having time-dependent covariates. The present study examines the nonlinear effects of the evolution of the covariates over time using penalized smoothing splines.

Cox’s proportional hazards model postulates that the hazard at time 𝑡 is the product of two components: ℎ(𝑡)=ℎ₀(𝑡)exp^𝐼_𝑖=1𝑏_𝑖𝑥_𝑖, (1) where 𝐛=(𝑏₁,…,𝑏_𝐼) is a vector of coefficients. The proportional hazards assumption is that the baseline hazard ℎ₀(𝑡) is a function of 𝑡 but does not involve the values of covariates. Several prognostic models for heart transplant survival data have been developed using Cox’s regression analysis, and the values of all covariates are determined at the time when the patient entered the study [7–9]. However, situations may exist in which the values of covariates change over the course of the study. A time-dependent model uses the follow-up data to estimate the effect of the evolution of the covariates over time. The relative hazard ℎ(𝑡)/ℎ₀(𝑡) then depends on time 𝑡, and thus the proportional hazards assumption is no longer satisfied [6, 10].

The time-dependent covariates 𝐗_𝑙^⟨𝑑⟩=(𝑥^⟨𝑑⟩_𝑙1,…,𝑥^⟨𝑑⟩_𝑙𝐼) are provided for patient no. 𝑑, where 𝑥^⟨𝑑⟩_𝑙1 is the midpoint of the 𝑙th time interval. Given the continuous survivor time, piecewise models arise from the partition of the time axis into disjointed intervals. Biganzoli et al. [11, 12] show that by treating the time interval as an input variable in a feed-forward neural network, it is possible to estimate smoothed discrete hazards as conditional probabilities of failure. To apply a generalized additive model (GAM), discretization of one-month or one-week intervals must be applied for the continuous survivor time with time-fixed covariates. However we cannot determine which discretization, one-month or one-week, should be applied; that is, the discretization is not initially unique. In the case of time-dependent covariates 𝐗_𝑙^⟨𝑑⟩, 𝑥^⟨𝑑⟩_𝑙1 is initially determined as the midpoint of the 𝑙th time interval for patient no. 𝑑. It is fairly straightforward to extend the model to survivor data with time-dependent covariates. Furthermore, by regarding a GAM as an extension of a partial logistic model (PLM), the unknown parameters can be estimated by maximizing the partial log-likelihood [13, 14].

We use the Stanford Heart Transplant data, which has been collected to model survival in patients. Although Cox’s proportional hazards model is not applicable in the case of time-dependent covariates, the survivor function can be estimated by taking ℎ₀(𝑡) in (1) to be the piecewise exponential hazard. Crowley and Hu [7], Aitkin et al. [8], and Lawless [9] used piecewise exponential models and plotted the survival function. Lagakos [15] also examined a graphical technique for assessing covariates in Cox’s proportional hazards model based on a permutation of the observed rank statistic. Most previous studies compared the hazard functions to assess the effect of transplantation on survival by fitting pretransplant and posttransplant data separately.

The difficulty is that there is no easily used measure of the difference between the transplanted and nontransplanted groups. Inferences must be based on a comparison of the estimated function. As Aitkin et al. [8] pointed out, there are always dangers in making inferences about the effect of treatment without adequate control groups. We thus provide an analysis that includes pretransplant and posttransplant data simultaneously as time-dependent covariates. It should be emphasized that patients who are not transplanted constitute a control group relative to patients who have undergone heart transplantation by the same covariates.

We use the 1977 version of the data, as given in Crowley and Hu [7], which is for 103 patients. As four of the transplanted patients have incomplete data on the mismatch score, our analysis is based on 99 patients to assess for what values of these covariates, if any, transplantation is likely to prolong survival. More than 30 percent of cases are censored. In these data, survival times are the number of days until death following a heart transplant, as in Lagakos [15]. A distinctive feature of the present problem is that some of the covariates are time-dependent (and possibly random). For example, Table 1 shows the values of covariates for transplant status (i.e., waiting time), age at transplant (in years), mismatch score (as time-dependent covariates), and previous open-heart surgery for patient no. 18. The previous surgery status does not change with time. In order to extend this setting, the covariate for transplant status is taken as an indicator (coded as 0 before the point of transplant and 1 after transplant). All the other time-dependent covariates are treated as being zero before transplant but changing from zero to the actual value of the particular covariate at the time of transplant. Patient no. 18 generated six observations. The proposed methods allow for simultaneous investigation of several covariates and provide estimates of the survival function as well as the significance.

Table 1 

Covariate values for patient no. 18.

2. Generalized Additive Models

By extending the PLM for the grouped data based on partial likelihood as introduced by Cox [16] and Efron [17], a PLM can be proposed for ungrouped data [13, 14] having time-dependent covariates for the discrete hazard rate ℎ_𝑙^⟨𝑑⟩ of patient no. 𝑑 at the time interval 𝑙: ℎPLM:ln^{𝑙⟨𝑑⟩}1−ℎ_𝑙^⟨𝑑⟩=𝛽₀+𝛽₁𝑥^⟨𝑑⟩_𝑙1+𝛽₂𝑥^⟨𝑑⟩_𝑙2+𝛽₃𝑥^⟨𝑑⟩_𝑙3+⋯+𝛽_𝐼𝑥^⟨𝑑⟩_𝑙𝐼. (2) In recent years, a variety of powerful techniques have been developed for exploring the functional form of effects. Here, GAM with smoothing splines proposed by Hastie et al. [18, 19] will be used by extending the generalized linear model (GLM) in McCullagh and Nelder [20], where the linear predictor in (2) is specified as a sum of smooth functions 𝑠(𝑥) with twice continuously differentiable functions of some or all of the covariates: ℎGAM:ln^{𝑙⟨𝑑⟩}1−ℎ_𝑙^⟨𝑑⟩=𝛽₀+𝑠₁𝑥^⟨𝑑⟩_𝑙1+𝑠₂𝑥^⟨𝑑⟩_𝑙2+𝑠₃𝑥^⟨𝑑⟩_𝑙3+⋯+𝑠_𝐼𝑥^⟨𝑑⟩_𝑙𝐼. (3)

The smooth functions in (3) can be represented as 𝑠₀(𝑥)=^𝑞0_𝑗=1𝛽_𝑗𝑏_0𝑗(𝑡),𝑠₁(𝑥)=^𝑞1_𝑗=1𝛽_𝑞0+𝑗𝑏_1𝑗(𝑥),⋮𝑠_𝐼(𝑥)=^𝑞𝐼_𝑗=1𝛽_{𝑞𝐼−1+𝑗}𝑏_𝐼𝑗(𝑥), (4) where 𝑞₁,𝑞₂,…,𝑞_𝐼 are the numbers of knots, and 𝜷=𝛽₀,𝛽₁,…,𝛽_𝑞0,𝛽_𝑞0+1,𝛽_𝑞0+2,…,𝛽_𝑞0+𝑞1,…,𝛽_{𝑞0+𝑞1+⋯+𝑞𝐼−1+1},𝛽_{𝑞0+𝑞1+⋯+𝑞𝐼−1+2},…,𝛽_{𝑞0+𝑞1+⋯+𝑞𝐼}. (5)

For time interval 𝑙 of patient no. 𝑑, we have the following definitions: 𝛿_𝑙^⟨𝑑⟩=𝛿1:patientno.𝑑diedattimeintervalforthe𝑙thclinicvisit,0:otherwise,_{′𝑙⟨𝑑⟩}=𝐯1:patientno.𝑑wascensoredattimeinterval𝑙,0:o.w,_{𝑙⟨𝑑⟩}=𝛿₁^⟨𝑑⟩,𝛿^′₁^⟨𝑑⟩,𝛿₂^⟨𝑑⟩,𝛿^′₂^⟨𝑑⟩,…,𝛿^⟨𝑑⟩_𝑙−1,𝛿^{′⟨𝑑⟩}_𝑙−1=(),𝐯0,0,…,0_{′𝑙⟨𝑑⟩}=𝐯_𝑙^⟨𝑑⟩,𝛿_𝑙^⟨𝑑⟩=0,0,…,0,𝛿_𝑙^⟨𝑑⟩, (6) where 𝐯_𝑙^⟨𝑑⟩ is the history of defaults and is censored for the first 𝑙−1 time intervals of patient no. 𝑑 and 𝐯^′_𝑙^⟨𝑑⟩=(𝐯_𝑙^⟨𝑑⟩,𝛿_𝑙^⟨𝑑⟩) is the same history extended to include 𝛿_𝑙^⟨𝑑⟩. Using the above model and notation, Tsujitani and Sakon [13] derived the full log-likelihood for all patients ln𝐿=ln𝐿(𝜷)+^𝑛𝑙_𝑑=1^𝑑_𝑙=1ln𝑝𝛿^′_𝑙^⟨𝑑⟩∣𝐯^′_𝑙^⟨𝑑⟩ (7) with the partial log-likelihood ln𝐿(𝜷)=^𝑛_𝑑=1^𝑙𝑑−1_𝑙=1ln1−ℎ_𝑙^⟨𝑑⟩+𝛿_𝑙^⟨𝑑⟩_𝑑lnℎ_𝑙^⟨𝑑⟩_𝑑+1−𝛿_𝑙^⟨𝑑⟩_𝑑ln1−ℎ_𝑙^⟨𝑑⟩_𝑑. (8) Although ln𝐿(𝜷) is not a log-likelihood in the usual sense, it possesses the usual asymptotic properties under fairly broad conditions, as proven in Andelsen and Gill [21]. To avoid overfitting, such models are estimated by penalized maximum likelihood ln𝐿(𝜷)=^𝑛_𝑑=1^𝑙𝑑−1_𝑙=1ln1−ℎ_𝑙^⟨𝑑⟩+𝛿_𝑙^⟨𝑑⟩_𝑑lnℎ_𝑙^⟨𝑑⟩_𝑑+1−𝛿_𝑙^⟨𝑑⟩_𝑑ln1−ℎ_𝑙^⟨𝑑⟩_𝑑+12^𝐼_𝑖=1𝜆_𝑖𝑠^′′_𝑖(𝑡)²𝑑𝑡, (9) where 𝜆_𝑖 are smoothing parameters that control the trade-off between the fit and the smoothness. The functions 𝑠_𝑖(𝑥) in (9) are represented by the 𝐵-spline basis functions 𝑏_𝑖(𝑥); see, for details, Tsujitani et al. [14].

Two model-fitting issues remain. The first concerns the selection of smoothing parameter 𝜆_𝑖 in (9). The optimum smoothing parameter choice is outweighed by the easy identification of a covariate’s functional form as well as the applicability of established inferential methods to short-term survival prediction. In order to select the smoothing parameters, the algorithm developed by Wood [22–24] can be applied by minimizing generalized cross validation (GCV) as an approximation to leave-one-out CV [23]. It should be noted that the leaving-one-out CV is to allow the deletion of only one observation. On the other hand, the ordinal 𝑣-fold CV divides the data randomly into 𝑣 groups so that their sizes are as nearly equal as possible. This partition should be made to avoid possible biases, as described in Zhang [25]. In many problems, the ordinal 𝑣-fold CV is, thus, unsatisfactory in several respects for time-dependent covariates. Applying this kind of data structure to the CV algorithm, we obtain insights into how the partition of data should be carried out. A natural extension of the 𝑣-fold CV algorithm by setting 𝑣=𝑛 is to allow the deletion of the patient with several observations; see, for details, Tsujitani et al. [14].

A second issue is the goodness-of-fit test of the model. After choosing the optimum smoothing parameters via the variant 𝑣-fold CV algorithm, the deviance allows us to test the goodness-of-fit: Dev=2ln𝐿_max−ln𝐿_𝑐, (10) where ln𝐿_𝑐 denotes the maximized partial log-likelihood under some current GAM and the log-likelihood for the maximum (full) model ln𝐿_max is zero. The deviance (10) is, however, not even approximately an 𝜒² distribution for the case in which ungrouped binary responses are available; see, for example, Collett [26], Landwehr et al. [27], and Tsujitani and Sakon [13]. The number of degrees of freedom required for the test for significance using the assumed 𝜒² distribution for the deviance is a contentious issue. No adequate distribution theory yet exists for the deviance. The reason for this is somewhat technical; for details, see Section 3.8 in Collett [2]. Consequently, the deviance on fitting a model to binary response data cannot be used as a summary measure of the goodness-of-fit of the model. Thus, bootstrapping is applied to the deviance (10) in order to obtain the goodness-of-fit; see, for details, Efron and Tibshirani [28] and Tsujitani et al. [14].

3. Example

As an initial model for the Stanford Heart Transplant data, we employ ModelI:𝑠𝑥₁+𝑥₂+𝑠𝑥₃+𝑠𝑥₄+𝑥₅. (11) GCV is only an approximation of leaving-one-out CV. Alternatively the variant 𝑣-fold CV is leaving-one-out CV based on each of 𝑛=99 patients to allow the deletion of the patient with several observations. By using variant 𝑣-fold CV and GCV for the initial model, the optimum smoothing parameters for GAM are determined as shown in Table 2. By using a backward elimination procedure, we obtain ModelII:𝑠𝑥₁+𝑥₂+𝑠𝑥₃+𝑥₅. (12) The likelihood ratio (LR) statistic based on deviance can be computed to test the significance of spline effects (i.e., nonlinearity). For example, the spline effect of “Midpoint(𝑥₁)” can also be tested by using ModelIII:𝑥₁+𝑥₂+𝑠𝑥₃+𝑥₅. (13) By comparing ModelII with III, the reduction in the value of deviance is Δ=4.59 with 1.85𝑑.𝑓. This is significant at the 10% level. The spline effect for “Age(𝑥₃)” is not significant. We thus obtain the final optimum GAM ModelIV:𝑠𝑥₁+𝑥₂+𝑥₃+𝑥₅ (14) with a variant 𝑣-fold score of 654.754.

Table 2 

Optimum smoothing parameters.

Figure 1 shows a histogram of the bootstrapped Dev(𝑏) for the optimum model. The bootstrap estimate of the 95th percentile Dev^∗ is Dev^∗=685.65. The comparison to Dev=639.77 of (10) suggests that the model fits the data.

Figure 1 

Histogram of the bootstrapped Dev(𝑏) for 𝐵=400.

Figure 2 shows the estimated contribution 𝑠(̂𝑥₁) of “Midpoint(𝑥₁)” to ln{ℎ_𝑙^⟨𝑑⟩/(1−ℎ_𝑙^⟨𝑑⟩)}, together with the ±2 standard deviation (SD) curves for the final optimum Model IV. The spline effects of 𝑥₁ are visualized in Figure 2. Figure 2 nicely shows that the spline function 𝑠(𝑥₁) of dying decreases initially as the midpoint 𝑥₁ increases. Subsequently, however, 𝑠(𝑥₁) is stably maintained after 500 midpoint. For the purpose of comparison, Figure 3 shows the estimated contribution 𝑠(̂𝑥₁) for GCV. From Figure 3, it is clear that the estimated 𝑠(̂𝑥₁) of 𝑥₁ is flat until 1500 and then tumbles because of too small smoothing parameter (i.e., overfitting), as shown in Table 2. So variant 𝑣-fold CV is superior to GCV. The analyses in this example are carried out using library{𝑚𝑔𝑐𝑣} in R.

Figure 2 

Estimated contribution 𝑠(̂𝑥1) of 𝑥1 (solid curve) and 𝑠(̂𝑥1)±2SD[𝑠(̂𝑥1)] (dashed curves).

Figure 3 

Estimated contribution 𝑠(̂𝑥1) of 𝑥1 (solid curve) and 𝑠(̂𝑥1)±2SD[𝑠(̂𝑥1)] (dashed curves) for GCV.

The survival function for our discretized situation is Pr(𝑙)=_{1≤𝑖≤𝑙}1−ℎ_𝑖^⟨𝑑⟩. (15) The average probability of survival at time interval 𝑙 for patient no. 𝑑 in group 𝑔 can be estimated as 𝑆_𝑔(𝑙)=1𝑛_𝑙^{𝑛[𝑔]𝑙[𝑔]}_𝑑=1Pr_𝑑^[𝑔](𝑙),𝑔=1,2,𝑙=1,2,…, (16) where 𝑛_𝑙^[𝑔] is the total number of patients at time interval 𝑙 in group 𝑔 and Pr_𝑑^[𝑔](𝑙) is the survival function Pr(𝑙) of patient no. 𝑑 at time interval 𝑙 in group 𝑔; see, for example, Thomsen et al. [29].

The data are analyzed to discover which values of the covariates are likely to be of benefit. We compare the results obtained using smoothing spline, partial logistic, Cox’s proportional hazards, and piecewise exponential models [7, 8]. The results of fitting the various models are summarized in Table 3. It is clear from Table 3 that(i)all covariates for the smoothing spline model are strongly significant (in particular, Crowley and Hu [7] suggested a quadratic effects of age) and(ii)there is little difference between Cox’s proportional hazard model and the piecewise exponential model. It should be noted that binary covariates in the model remain linear.

Table 3 

𝑃 values for the significance test of covariates.

As shown in Aitkin et al. [8, Figure 2], it is more appropriate to compare survivorship functions if the hazards are not proportional. One point of interest is a comparison of survival experience of transplanted and nontransplanted patients. Our proposal for comparing the survival function is to use the estimated survival function for only 41 heart transplanted patients who died to assess the efficacy of transplantation and the effects of covariates by means of modeling the change in hazard at transplantation by using (15) and (16). Our particular interest is the effect of waiting time on posttransplant survival according to several models. In Figure 4, two time periods are used (group 1: up to 20 days; group 2: longer than 20 days). Figure 4 shows a comparison of the estimated survival function. The estimated survival functions based on the smoothing spline suggest that patients with a short waiting time face a greater early risk than those who had a longer waiting time. However the estimated survival functions based on piecewise exponential models cannot reveal the difference between a short and long waiting times. Our method provides an alternative to Arjas’ [10] suggestion of comparing separate estimates of cumulative hazard based on the levels of the waiting time. Although Arjas [10] did not include waiting time as a covariate in Cox’s proportional hazard model because of nonproportionality issues, we used transplant status (i.e., waiting time), which is strongly significant for the smoothing spline model according to the results shown in Table 3.

Figure 4 

Survival function describing the effect of the waiting time for 41 heart transplanted patients who died.

A fundamentally different type of analysis was suggested by Crowley and Hu [7] to investigate the effect of transplantation with a low mismatch score. They pointed out that transplantation may be beneficial for younger patients only based on regression coefficients for Cox’s proportional hazards model, but our conclusion can be derived by graphical analysis as well as significance testing of covariates. Defining a low mismatch score as less than or equal to one for all 29 heart transplanted patients [7], Figure 5 shows a graphical comparison of the estimated survival function for two groups, namely, the younger patients (less than 50 years old at acceptance) and older patients (greater than or equal to 50 at acceptance). From Figure 5, it is clear that older patients face a greater early risk than younger patients; see, for details, Crowley and Hu [7, Chapter 5] with respect to the cutpoints for low mismatch score as less than or equal to one and the younger patients as less than 50 years old. Kalbfleish and Prentice [30, Section  4.6.3] estimated the cutpoint for age, based on all 65 transplanted patients, as 46.2. Figure 6 shows a graphical comparison of the estimated survival function for two groups, namely, the younger patients (less than or equal 46 years old at acceptance) and older patients (greater than 46 at acceptance). As Kalbfleish and Prentice point out, transplantation is beneficial for younger patients.

Figure 5 

Survival function describing the effect of the age of transplantation for patients with a low mismatch score.

Figure 6 

Survival function describing the effect of the age for all 65 transplanted patients.

4. Conclusion

We confined our attention to time-dependent covariates. Allowing covariates to vary over the duration of the study not only enabled us to study time-varying risk factors, but also provided a flexible way for modeling censored survival data using penalized smoothing splines. We illustrated the procedures using data of the Stanford Heart Transplant data.

By introducing the maximum likelihood principle into GAM,(i)we could visualize the spline effects of the midpoint of the time interval;(ii)the smoothing parameters could be selected by using variant 𝑣-fold CV;(iii)the goodness-of-fit of GAM could be tested based on bootstrapping;(iv)the estimated average probabilities of survival enabled us to investigate the effect of transplantation with a low mismatch score for two groups, namely, the younger and older patients.