Survival Data Analysis with Time-Dependent Covariates Using Generalized Additive Models

Abstract

We discuss a flexible method for modeling survival data using penalized smoothing splines when the values of covariates change for the duration of the study. The Cox proportional hazards model has been widely used for the analysis of treatment and prognostic effects with censored survival data. However, a number of theoretical problems with respect to the baseline survival function remain unsolved. We use the generalized additive models (GAMs) with B splines to estimate the survival function and select the optimum smoothing parameters based on a variant multifold cross-validation (CV) method. The methods are compared with the generalized cross-validation (GCV) method using data from a long-term study of patients with primary biliary cirrhosis (PBC).

1. Introduction

Several prognostic models for PBC data have been developed using the Cox proportional hazards model, and the values of all covariates were determined at the time when the patient entered the study [1]. However, situations may exist in which the values of covariates change for the duration of the study. The time-dependent model uses follow-up data to estimate the effect of the evolution of the covariates during the course of the disease; see, for example, Cox [2], Altman and Stavola [3] and Collett [4].

Let 𝑡 be a continuous lifetime variable and 𝐱=(𝑥₁,…,𝑥_𝐼) a vector of time-fixed covariates. The Cox’s proportional hazards model postulates that the hazard at time 𝑡 is the product of two components [5, 6]ℎ(𝑡;𝐱)=ℎ₀(𝑡)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 𝐱 which are measured at the beginning of an interval to predict short-term survival.

We investigate PBC data for 312 patients who were seen at the Mayo Clinic and were monitored for the duration of the study, as described in Murtaugh et al. [7] and Therneau and Grambsch [8]. The Cox proportional hazards model was developed based on the relationship between survival and the patient characteristics observed when the patient entered the study. The precision of time-fixed models used in PBC is rather low, partly because these models are based on data for which the covariates were measured at the time when the patient entered the study.

For the analysis of data with time-dependent covariates, however, the survivor function for any individual depends on time 𝑡 and the baseline hazard function. This means that the survivor function cannot be expressed as a power of the baseline survivor function and is generally difficult to obtain for any individual; see, for example, Kalbfleish and Prentice [9] and Marubini and Valsecchi [10]. The Mayo updated model (e.g., [7]), and the European new version model (e.g., [3, 11, 12]) have been commonly used to improve the accuracy of survival predictions as a function of covariates measured at any time during the course of the disease. In the present article, we propose the variant multifold CV method for GAM when choosing the optimum smoothing parameters in order to estimate the survival function and predict the short-term survival (say, for the following six months) at any time during the course of the disease.

Another useful idea in our analysis is the concept of competing risk. There is “liver transplantation” in PBC data as competing risk. Competing risk has been treated as censored data. By adding the liver transplantation as one of time-dependent covariate, one can test the significance of liver transplantation.

2. Model Building

By extending the Cox proportional hazard model (1), a flexible survival model has been examinedℎ(𝑡;𝐱)=ℎ₀(𝑡)exp^𝐼_𝑖=1𝑠_𝑖𝐱_𝑖, (2) where 𝑠_𝑖(𝐱_𝑖) is a spline function for the covariate 𝐱_𝑖 [13–15]. The proportional hazard model (2) used the time-fixed values of covariates as shown in Dickson et al. [1]. The estimates of hazard ratio by relative survival regression model [16] with time-dependent covariates are compared with that of Cox proportional hazard model. A new approach [17, 18] is proposed with PBC data, aiming to capture nonlinear patterns of bilirubin time courses and their relationship with survival time of patients. However, because most patients with PBC make repeated visits to the clinic, it is natural to ask the optimum timing of liver transplantation by predicting short-term survival at any time in the course of the disease.

The time-dependent covariates 𝐗_𝑙^⟨𝑑⟩=(𝑎_𝑙^⟨𝑑⟩,𝑥^⟨𝑑⟩_𝑙1,…,𝑥^⟨𝑑⟩_𝑙𝐼) are provided for patient #𝑑, where 𝑎_𝑙^⟨𝑑⟩ is the midpoint of time interval [𝑡_𝑙^⟨𝑑⟩,𝑡^⟨𝑑⟩_𝑙+1] for the 𝑙th clinic visit. The event times may be subject to the usual random censoring. Then only the minimum of survival and censoring time with censoring indicator𝛿_𝑙^⟨𝑑⟩=1∶patient#𝑑diedattime-intervalforthe𝑙thclinicvisit0∶otherwise (3)

are observed. The relative hazard ℎ(𝑡)/ℎ₀(𝑡) then depends on time 𝑡, and thus the proportional hazards assumption is no longer satisfied, as described in Altman and Stavola [3] and Arjas [19].

For example, Table 1 shows the values of age, prothrombin time, and bilirubin as time-dependent covariates for dead patient #9; for details, see Table 4 in Murtaugh et al. [7]. Patients were scheduled to return for further observations at six months, 12 months, and yearly. Thus, 𝑛=312 patients generate 1945 observations in total. The covariates values for each patient were allowed to vary with the time interval for the 𝑙th clinic visit.

Table 1 

Values of covariates for deceased patient #9.

A grouped version of Cox’s proportional hazard model with time-fixed covariates has been considered in the framework of discrete grouped data for the feed-forward neural network. Given the continuous survivor time, piecewise models arise from the partition of the time axis into disjointed intervals. Biganzoli et al. [20, 21] 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. Biganzoli et al. [20] also pointed out that an advantage of this kind of data structure is the possibility of straightforward use of time-dependent covariates since each subject is represented, for each observation interval, by one input vector which can change across intervals. In order to apply this neural network approach, which is called partial logistic regression models [20], discretization of one-month or one-week intervals must be applied for the continuous survivor time with time-fixed covariates. We cannot determine which discretization, one-month or one-week intervals, must be applied; that, is the discretization is not originally unique. For the data in Table 1, however, the choice of discretization of the time axis for the partial logistic regression model is generally determined by clinical relevance, possibly according to the scheduled time intervals between follow-up visits.

The primary goal of the present study is to predict short-term survival in patients on the basis of measurements of several characteristics having time-dependent covariates 𝐗_𝑙^⟨𝑑⟩ for the purpose of facilitating the decision as to when to undertake liver transplantation. Based on partial logistic model due to Cox [22] and Efron [23] for the grouped data, Tsujitani and Sakon [24] have proposed a partial logistic model with a discrete hazard rate ℎ_𝑙^<𝑑> for ungrouped data having time-dependent covariatesℎln^{𝑙⟨𝑑⟩}1−ℎ_𝑙^⟨𝑑⟩=𝛽₀+𝛾𝑡_𝑙^⟨𝑑⟩+𝛽₁𝑥^⟨𝑑⟩_𝑙1+𝛽₂𝑥^⟨𝑑⟩_𝑙2+⋯+𝛽_𝐼𝑥^⟨𝑑⟩_𝑙𝐼,𝑑=1,2,….,𝑛, (4) where 𝜷=(𝛽₁,…,𝛽_𝐼) is a vector of coefficients. The modeled response is the logit of hazard rate, and the logit is linear in the covariates. However, this assumption is violated when covariate effects are best represented by smooth, nonlinear function. In recent years, a variety of powerful techniques have been developed for exploring the function form of effects. We examine here a flexible survival model GAM that does not require linearity of the covariate function by extending a generalized linear model (GLM); see, for example, Hastie and Tibshirani [13] and McCullagh and Nelder [25]. By identification of nonlinear covariate effects, we can estimate more accurately a patient’s prognosis and thus determine a liver transplant based on prediction of short-term survival.

The linear predictor in (4) is specified as a sum of smooth functions 𝑠(𝐱) with twice continuous derivatives of some or all of the covariates for the discrete hazard rate ℎ_𝑙^<𝑑> of patient #𝑑 at the time interval 𝑙ℎln^{𝑙⟨𝑑⟩}1−ℎ_𝑙^⟨𝑑⟩=𝛽₀+𝑠₀𝑡_𝑙^⟨𝑑⟩+𝑠₁𝑥^⟨𝑑⟩_𝑙1+𝑠₂𝑥^⟨𝑑⟩_𝑙2+⋯+𝑠_𝐼𝑥^⟨𝑑⟩_𝑙𝐼. (5)

The smooth functions in (5) can be represented as𝑠₀(𝑥)=^𝑞0_𝑗=1𝛽_𝑗𝑏_0𝑗(𝑡),𝑠₁(𝑥)=^𝑞1_𝑗=1𝛽_𝑞0+𝑗𝑏_1𝑗(𝑥),…𝑠_𝐼(𝑥)=^𝑞𝐼_𝑗=1𝛽_{𝑞𝐼−1+𝑗}𝑏_𝐼𝑗(𝑥), (6)where 𝑞₁,𝑞₂,…,𝑞_𝐼 are the numbers of knots, and𝜷=𝛽₀,𝛽₁,…,𝛽_𝑞0,𝛽_𝑞0+1,𝛽_𝑞0+2,…,𝛽_𝑞0+𝑞1,…,𝛽_{𝑞0+𝑞1+⋯+𝑞𝐼−1+1},𝛽_{𝑞0+𝑞1+⋯+𝑞𝐼−1+2},…,𝛽_{𝑞0+𝑞1+⋯+𝑞𝐼}, (7)𝐙=1𝑏₀₁𝑥₁₁…𝑏_0𝑞1𝑥₁₁⋅𝑏_{𝐼−1,𝑞0+𝑞1+⋯+𝑞𝐼−1+1}𝑥_𝐼−1,1…𝑏_{𝐼,𝑞0+𝑞1+⋯+𝑞𝐼}𝑥_𝐼11𝑏₀₁𝑥₁₂…𝑏_0𝑞1𝑥₁₂⋅𝑏_{𝐼−1,𝑞0+𝑞1+⋯+𝑞𝐼−1+1}𝑥_𝐼−1,2…𝑏_{𝐼,𝑞0+𝑞1+⋯+𝑞𝐼}𝑥_𝐼2⋅⋅…⋅⋅……⋅1𝑏₀₁𝑥₁₁…𝑏_0𝑞1𝑥_1𝑛⋅𝑏_{𝐼−1,𝑞0+𝑞1+⋯+𝑞𝐼−1+1}𝑥_{𝐼−1,𝑛}…𝑏_{𝐼,𝑞0+𝑞1+⋯+𝑞𝐼}𝑥_𝐼𝑛. (8)The functions 𝑏_𝑖𝑗(𝑥_𝑖𝑗) in the matrix 𝐙 are B-spline basis functions. Thus, (5) can be rewritten asℎln^{𝑙⟨𝑑⟩}1−ℎ_𝑙^⟨𝑑⟩=𝐙𝜷_{𝑛×𝑞𝑞×1};𝑞=^𝐼_𝑖=0𝑞_𝑖, (9) where B-splines with 10 interiors knots will be used for each continuous covariate. The number of knots is arbitrary but appears to have little effect on the results, provided that the number is not too small, as described in Gray [15].

At the time interval for the 𝑙th clinic visit of patient #𝑑, we define𝛿_𝑙^{′⟨𝑑⟩}=𝐯1∶patient#𝑑wascensoredattime-intervalforthe𝑙thclinicvisit0∶otherwise,_{𝑙⟨𝑑⟩}=𝛿₁^⟨𝑑⟩,𝛿₁^{′⟨𝑑⟩},𝛿₂^⟨𝑑⟩,𝛿₂^{′⟨𝑑⟩},…,𝛿^⟨𝑑⟩_𝑙−1,𝛿^{′⟨𝑑⟩}_𝑙−1=(),𝐯0,0,…,0_{𝑙′⟨𝑑⟩}=𝐯_𝑙^⟨𝑑⟩,𝛿_𝑙^⟨𝑑⟩=0,0,…,0,𝛿_𝑙^⟨𝑑⟩, (10)

where 𝐯_𝑙^⟨𝑑⟩ is the history of dead and censored of time intervals for the first 𝑙th clinic visit of patient #𝑑, and 𝐯_𝑙^{′⟨𝑑⟩}=(𝐯_𝑙^⟨𝑑⟩,𝛿_𝑙^⟨𝑑⟩) is the same history extended to include 𝛿_𝑙^<𝑑>. Tsujitani and Sakon [24] derived the full log likelihood for all patientsln𝐿=ln𝐿(𝜷)+^𝑛𝑙_𝑑=1^𝑑_𝑙=1ln𝑝𝛿_𝑙^{′⟨𝑑⟩}∣𝐯_𝑙^{′⟨𝑑⟩} (11) with partial log likelihoodln𝐿(𝜷)=^𝑛_𝑑=1^𝑙𝑑−1_𝑙=1ln1−ℎ_𝑙^⟨𝑑⟩+𝛿_𝑙^⟨𝑑⟩_𝑑lnℎ_𝑙^⟨𝑑⟩_𝑑+1−𝛿_𝑙^⟨𝑑⟩_𝑑ln1−ℎ_𝑙^⟨𝑑⟩_𝑑. (12) The unknown parameters 𝜷 in (9) can thus be estimated by maximizing the partial log likelihood (12), which is the log likelihood for the independent Bernoulli trial. Although ln𝐿(𝜷) is not a log likelihood in the usual sense, it possesses the usual asymptotic properties under fairly broad conditions; see, for example, Andelsen and Gill [26].

To avoid overfitting, such models are estimated by penalized maximum likelihoodln𝐿(𝜷)=^𝑛_𝑑=1^𝑙𝑑−1_𝑙=1ln1−ℎ_𝑙^⟨𝑑⟩+𝛿_𝑙^⟨𝑑⟩_𝑑lnℎ_𝑙^⟨𝑑⟩_𝑑+1−𝛿_𝑙^⟨𝑑⟩_𝑑ln1−ℎ_𝑙^⟨𝑑⟩_𝑑+12^𝐼_𝑖=1𝜆_𝑖𝑠_𝑖^′′(𝑡)²𝑑𝑡, (13) where 𝜆_𝑖 are smoothing parameters that control the tradeoff between the fit and the smoothness, and 𝑠_𝑖^′′ is the twice derivative of 𝑠_𝑖 with respect to 𝑡. The advantage of penalized estimates is enlightened in Wood ([27], Section  4.1).

Two model-fitting issues remain. The first concerns the selection of smoothing parameter 𝜆_𝑖 in (13). The careful smoothing parameter choice is outweighed by the easy identification of a covariate’s functional form, and the applicability of established inferential methods to short-term survival prediction. In order to select the smoothing parameters, the algorithm due to Wood [27–29] can be used by minimizing GCV as an approximation to leaving-one-out CV. For example, however, the dead patient #9 generated seven observations as shown in Table 1. Patients were scheduled to return for further observations at six months, 12 months, and yearly. It should be noted that this patient generated seven observations. Thus, 𝑛=312 patients generate 1945 observations in total.

We propose a natural extension of 𝑣-fold CV algorithm by “leaving-one-out” CV based on each 𝑛=312 patients. The ordinal 𝑣-fold CV divides the data randomly in 𝑣 groups so that their sizes are as nearly equal as possible. The partition should be made to avoid possible biases, as described in Zhang [30]. 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 done. A natural extension of 𝑣-fold CV algorithm by setting 𝑣=𝑛 is to allow the deletion of the patient with several observations. The variant 𝑣-fold CV is given as follows:

Step 1. Split the original sample 𝐗={𝐗^⟨1⟩,…,𝐗^⟨𝑑⟩,…,𝐗^⟨𝑛⟩} into 𝑛 parts 𝐗={𝐗^⟨1⟩|…|𝐗^⟨𝑑⟩|…|𝐗^⟨𝑛⟩}, where 𝐗^⟨𝑑⟩={𝐗₁^⟨𝑑⟩,…,𝐗_𝑙^⟨𝑑⟩_𝑑}.

Step 2. Fit the model to 𝐗^∗={𝐗^⟨1⟩|…|𝐗^{⟨𝑑−1⟩}|𝐗^⟨𝑑+1⟩|…|𝐗^⟨𝑛⟩} with the 𝑑-th subject (i.e., patient) deleted of the data, and predict ∑^𝑙𝑑−1_𝑙=1ln(1−ℎ_𝑙^⟨𝑑⟩)+𝛿_𝑙^⟨𝑑⟩_𝑑lnℎ_𝑙^⟨𝑑⟩_𝑑+(1−𝛿_𝑙^⟨𝑑⟩_𝑑)ln(1−ℎ_𝑙^⟨𝑑⟩_𝑑) for the deleted 𝑑th sample 𝐗^⟨𝑑⟩.

Step 3. Do the above for 𝑑=1,2,…,𝑛 and combine the CV estimates CV=^𝑛_𝑑=1^𝑙𝑑−1_𝑙=1ln1−ℎ_𝑙^⟨𝑑⟩+𝛿_𝑙^⟨𝑑⟩_𝑑lnℎ_𝑙^⟨𝑑⟩_𝑑+1−𝛿_𝑙^⟨𝑑⟩_𝑑ln1−ℎ_𝑙^⟨𝑑⟩_𝑑. (14)

A second issue is the goodness-of-fit test of the model. After choosing the optimum smoothing parameters via 𝑣-fold CV algorithm, the deviance allows us to test the goodness of fit Dev=2ln𝐿_max−ln𝐿_𝑐, (15) 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 given by (15) is, however, not even approximately a 𝜒² distribution for the case in which ungrouped binary responses are available; see, for example, Landwehr et al. [31] and Tsujitani and Koshimizu [32] and Collett [33]. 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 exists for the deviance. The reason for this is somewhat technical; for details, see Section  3.8 in Collett [33]. Consequently, the deviance on fitting a model to binary response data cannot be used as a summary measure of the goodness-of-fit test of the model.

Based on the above discussion, we employ bootstrapping to the deviance of (15) in order to obtain the goodness-of-fit test due to Efron and Tibshirani [34].

Step 1. Generate 𝐵 bootstrap samples 𝐗^∗={𝐗_𝑙^⟨1⟩∗,…,𝐗_𝑙^{⟨𝑛⟩∗}} from the original sample 𝐗. Let 𝐗^∗(𝑏) denote the 𝑏th bootstrap sample.

Step 2. For the bootstrap sample 𝐗^∗(𝑏), the deviance of (15) is computed by Dev(𝑏)=2[ln𝐿_𝑓−ln(𝐗^∗(𝑏);𝜷^∗(𝑏))].
This process is repeated independently 𝐵 times, and the computed values are arranged in ascending order.

Step 3. The value of the 𝑗th order statistic Dev^∗ of the 𝐵 replications can be taken as an estimate of the quantile of order 𝑗/(𝐵+1).

Step 4. The estimate of the 100(1−𝛼)-th percentile (i.e., 𝛼% critical point) of Dev^∗ is used to test the goodness-of-fit of a model having a specified significance level 𝛼=1−𝑗/(𝐵+1). The value of deviance of (15) being greater than the estimate of the percentile indicates that the model fits poorly.

3. Results

The survival function for our discretized situation is𝑆𝑡_𝑙=_{1≤𝑖≤𝑙}1−ℎ_𝑖^⟨𝑑⟩. (16) The conditional probability Pr(𝑡,𝑡+Δ𝑡) of survival over a short-time interval Δ𝑡 (say, six months) after time 𝑡 can be estimated as()=𝑆()Pr𝑡,Δ𝑡𝑡+Δ𝑡𝑆(𝑡). (17)

By using variant 𝑣-fold CV and GCV, the optimum smoothing parameters for GAM are determined as shown in Table 2. Table 3 summarizes the 𝑃 values to test the nonparametric effects of covariates for the model s(time) + s(age) + s(pro) + s(bili) with the optimum smoothing parameters and GCV. From Table 3, all covariates are highly significant for GCV; however, time is not significant for variant 𝑣-fold CV. GCV is only the approximation of leaving-one-out CV. Furthermore, the variant 𝑣-fold CV is leaving one-out CV based on each 𝑛=312 patients. So variant 𝑣-fold CV is better than GCV. For the purpose of comparison, we included the results in the case using GCV in Tables 2 and 3. Table 4 shows the results of a number of models that were fit to the data.

Table 2 

Optimum smoothing parameters.

Table 3 

Test of significance for the covariates (𝑃 values).

Table 4 

Ranks and 𝐹 values for nonparametric effects.

The likelihood ratio (LR) statistic based on deviance can be conducted to test whether the spline effect provides a significantly better fit than a linear effect. Table 5 shows the test of significance for spline effects based on the models in Table 4. It is clear from Table 5 that the spline effects of prothrombin time and bilirubin are strongly significant. No spline provides a significantly better fit than a linear model for age. Thus, we accept the final model: age + s(pro) + s(bili) with the edf(effective degrees of freedom) 3.954 and 4.477 for s(pro) and s(bili), respectively. We used likelihood ratio test instead of information criteria as a valid alternative approach for model selection. From Table 4, however, it is found that the same final model is selected by using AIC.

Table 5 

Test of significance for spline effects.

Figure 1 shows the histogram of the bootstrapped Dev(𝑏) for the optimum model with 𝐵=400. The bootstrap estimate of the 95th percentile (i.e., 5% critical point) Dev^∗ is Dev^∗=734.59. Comparison to Dev=663.65 of (15) suggests that the model fits the data.

Figure 1 

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

Figure 2 shows the prediction of the probability of surviving beyond the next six months for dead patient #9. For the purpose of comparison, the results obtained using partial logistic regression, the Mayo updated, and the European new version models are also provided. Figure 2 also indicates that the six-month survival probability predicted by GAM are lower than those predicted by the other models. Because the patient #9 died, the lower predicted probabilities are better. The conditional probability Pr(𝑡,Δ𝑡) of survival over a short time interval Δ𝑡 (say, six months) after time 𝑡 during the course of the disease can be predicted from data collected for censored and dead data.

Figure 2 

Probability of survival over the next six months using the four models with respect to dead patient #9.

For the graphical representation, the individual probabilities for predicted survival are averaged in order to compare the Mayo updated model, the European new version model, the partial logistic regression model, and GAM. We can predict the probability of survival over the following six months using the four models with respect to data and censored data out of 312 patients. For the group (𝑔=1) and the censored group (𝑔=2), the probability of surviving over the next Δ𝑡 months is denoted by Pr_𝑑^[𝑔](𝑙,𝑙+Δ𝑡) for the 𝑙-th clinic visit of patient #𝑑. The average probability of survival over the next Δ𝑡 months for the 𝑙-th clinic visit of patient #𝑑 in group 𝑔 can be estimated as𝑆_𝑔(𝑙)=1𝑛_𝑙^{𝑛[𝑔]𝑙[𝑔]}_𝑑=1Pr_𝑑^[𝑔](𝑙,𝑙+Δ𝑡),𝑔=1,2,𝑙=1,2,…,𝐿, (18) where 𝑛_𝑙^[𝑔] is the total number of patients for the 𝑙-th clinic visit in group 𝑔, and Pr_𝑑^[𝑔](𝑙) is the survival function Pr(𝑙) of patient #𝑑 at the time interval 𝑙 in group 𝑔; see, for example, Markus et al. [35], Marubini and Valsecchi [10], and Thomsen et al. [36].

Figure 3 shows a comparison of the probability of survival over the next six months using the four models with respect to dead and censored data among all 312 patients. The figure clarifies that(i)for the case of dead data, the six-month survival probabilities predicted by GAM are lower than those predicted by the other models, and,(ii)for the case of censored data, the difference among the four models is very small.

Figure 3 

Probability of survival over the next six months using the four models with respect to dead and censored data among all 312 patients.

Figure 4 also shows the box and whisker plots of probability of survival over the next six months using GAM with respect to dead data among all 312 patients. It should be noted that the variance of probabilities of survival over the next six months is much higher in the fourth clinic visits than in other clinic visits. Another useful idea in our analysis is the concept of competing risk. There is “liver transplantation” in PBC data as competing risk. Competing risk has been treated as censored data. By adding 𝑥^⟨𝑑⟩_𝑙4 as one of time-dependent covariate for the liver transplantation, one can test the significance of liver transplantation. The covariate for liver transplantation is taken as a binary variable (codes 0 before liver transplantation, 1 at liver transplantation) as shown in Giorgi and Gouvernet [16] and Crowley [37]. Table 6 shows the three types for the combination of “censored” and “liver transplantation.” Table 7 shows the values of covariates for liver-transplanted patient #5.

Table 6 

Three types of “censored” and “liver transplantation.”

Table 7 

Values of covariates for dead patient #5.

Figure 4 

Box and whisker plots of probability of survival over the next six months using GAM with respect to dead data among all 312 patients.

In order to test the significance of “liver transplantation,” we consider two models:

Model IV: age + s(pro) + s(bili).

Model IV′: age + s(pro) + s(bili) + liver transplantation. The values of deviance and d.f. are given in Table 8. The reduction in the value of deviance is 10.657=663.65−653.00 on 0.973 d.f. This is significant at the 1% level.

Table 8 

Test of significance for covariates using model (19).

For the purpose of the comparison, the hazard of the cumulative incidence function (CIF) may be modeled in the presence of competing risks. The model is based on 𝛾(𝑡;𝐱)=𝛾₀(𝑡)exp^𝐼_𝑖=1𝛽_𝑖𝑥_𝑖, (19) where 𝑡 is the time of the last observation (not the midpoint at the time interval 𝑙), 𝛾 is the hazard of the subdistribution, and 𝛾₀ is the baseline hazard of the subdistribution ([38], Section  6.2). The 𝑃 values are summarized in Table 8 to test the significance for covariates using the model (19). From Tables 3 and 8, there is little difference between our method and the CIF.

4. Discussion

In this paper, we introduced the probabilistic interpretation of GAM and constructed the maximum likelihood principle of GAM for the analysis of survival data having time-dependent covariates. We proposed the information criterion based on the variant 𝑣-fold CV when choosing the optimal smoothing parameters in application of GAM. Introducing the maximum likelihood principle into GAM, the deviance allows us to test the goodness-of-fit of GAM. The proposed methods were illustrated by comparing the probability of survival over the next six months using the Mayo-updated model, the European new version model, the partial logistic regression model, and GAM with respect to dead and censored data among PBC data. We expect that flexible methods for modeling survival data with time-dependent covariates using machine learning theory such that support vector machine will be very useful in this real-world contexts; see, for example, Hastie et al. [39]. Furthermore, smoothing spline ANOVA models by Gu [40] will enable us to include the interactions between the covariates.

We assume that there is only one cause of failure; that is, the event is allowed to occur only once for each patient. However, there is increasing interest to apply survival data sets with multiple events per patient [8, 41]. Wei et al. [42] analyzed bladder cancer data by modeling marginal distributions of multivariate failure time with proportional hazards models. The model may violate the proportional hazards assumption, even when the overall data set does not (Table 9). By modifying such as𝛿_𝑙^⟨𝑑⟩=1:recurrenteventforpatient#𝑑occurredatthetime-interval𝑙0∶otherwise. (20)

Table 9 

Deviance and d.f.

Equation (11) becomesln𝐿(𝜷)=^𝑛𝑙_𝑑=1^𝑑_𝑙=1𝛿_𝑙^⟨𝑑⟩_𝑑lnℎ_𝑙^⟨𝑑⟩_𝑑+1−𝛿_𝑙^⟨𝑑⟩_𝑑ln1−ℎ_𝑙^⟨𝑑⟩_𝑑. (21)

Thus, the ideas presented in this paper can be extended to identification of prognostic factors relative to survival time in the case that the same event may recur during a follow-up study, and covariates values change with time.