Analyzing Competing Risk Data Using the R timereg Package

In this paper we describe ﬂexible competing risks regression models using the comp.risk() function available in the timereg package for R based on Scheike et al. (2008). Regression models are speciﬁed for the transition probabilities, that is the cumulative incidence in the competing risks setting. The model contains the Fine and Gray (1999) model as a special case. This can be used to do goodness-of-ﬁt test for the subdistribution hazards’ proportionality assumption (Scheike and Zhang 2008). The program can also construct conﬁdence bands for predicted cumulative incidence curves. We apply the methods to data on follicular cell lymphoma from Pintilie (2007), where the competing risks are disease relapse and death without relapse. There is important non-proportionality present in the data, and it is demonstrated how one can analyze these data using the ﬂexible regression models.


Introduction
Competing risks data often arise in biomedical research when subjects are at risk of failure from K different causes. When one event occurs, it precludes the occurrence of any other event. In cancer studies, one common example of competing risks involves disease relapse and death in remission. The cumulative incidence curve, i.e., the probability of failure of a specific type is a useful summary curve when analyzing competing risks data. Unfortunately this is not widely known in the biomedical world, and a very common error is that people report one minus the Kaplan-Meier estimate for each competing cause as a probability of cause-specific free survival. This is not a correct procedure and this estimator overestimates the incidence rates of a particular cause in the presence of all other competing causes (see Klein et al. 2001 for details).
The aim of this work is to estimate and model the cumulative incidence probability of a specific cause of failure. Estimating and modelling the cause-specific hazards has been considered as a standard approach for analyzing competing risks data. Assuming two types of failures k = 1, 2, the cumulative incidence function for cause 1 given a set of covariates x is given by where T is the failure time, indicates the cause of failure and λ k (t; x) is the hazard of the kth cause failure conditional on x, which is defined as Here, the cause-specific hazards for all causes need to be properly modeled. Cox's proportional hazards model is the most popular regression model in survival analysis and here the hazard function is given by where λ k0 (t) is a cause-specific baseline and β are regression coefficients. Using Cox's regression model to model the cause-specific hazards with the purpose of estimating the cumulative incidence function (1) was considered by Lunn and McNeil (1995) and Cheng et al. (1998). Shen and Cheng (1999) considered Lin and Ying's special additive model for the cause-specific hazards and Zhang (2002, 2003) considered a flexible Cox-Aalen model. The latter model allows some covariates to have time-varying effects. Modelling of the cause-specific hazards gives a complex nonlinear modelling relationship for the cumulative incidence curves. It is therefore hard to summarize the covariate effect and hard to identify the time-varying effect on the cumulative incidence function for a specific covariate. Recently, it has been suggested to directly model the cumulative incidence function. Fine and Gray (1999, FG) developed a direct Cox regression approach to model the subdistribution hazard function of a specific cause. The cumulative incidence function based on the FG model is given by where Λ 1 (t) is an unknown increasing function and β is a vector of regression coefficients. FG proposed using an inverse probability of censoring weighting technique to estimate β and Λ 1 (t). This approach is implemented in the crr() function in the cmprsk package (Gray 2010) for R (R Development Core Team 2010).
Recently, we considered a class of flexible models of the form where h and g are known link functions and α(t) and γ are unknown regression coefficients (see Scheike et al. 2008, SZG). FG's proportional regression model, Lin and Ying's special additive model and Aalen's full additive regression model are special sub-models of our model. Any link function can be considered and used here. In this study we focus on two classes of flexible models: proportional models and additive models The regression coefficients α(t) and γ are estimated by a simple direct binomial regression approach. We have developed a function, comp.risk(), available in the R package timereg, that implements this approach. In addition we have proposed a useful goodness-of-fit test to identify whether time-varying effect is present for a specific covariate.
In medical studies physicians often wish to estimate the predicted cumulative incidence probability for a given set of values of covariates. The predict() function of timereg computes the predicted cumulative incidence probability and an estimate of its variance at each fixed time point, and constructs (1 − α)100% simultaneous confidence bands over a given time interval.
One further advantage is that the software can deal with cluster structure, see Scheike et al. (2010).
The estimation procedure and goodness-of-fit test will be presented in Section 2. In Section 3 we will show how the comp.risk() function in the R package timereg can be used to fit our newly proposed flexible models (3) and (4)

Estimation
Let T i and C i be the event time and right censoring time for the ith individual, respectively.
i ∈ {1, . . . , K} indicates the cause of failure. LetT i = min(T i , C i ) and ∆ i = I(T i ≤ C i ). We observe n independent identically distributed (i.i.d.) realizations of {T i , ∆ i , ∆ i i , X i , Z i } for i = 1, . . . , n, where X i = (1, X i1 , . . . , X ip ) and Z i = (Z i1 , . . . , Z iq ) are associated covariates. We assume that (T i , i ) are independent of C i given covariates. Let N i (t) = I(T i ≤ t, i = 1) be the underlying counting processes associated with cause 1, which are not observable for all t. However, ∆ i N i (t) are computable for all t and we can show that where G(t|X, Z) is the survival distribution for the censoring time given the covariates. We therefore considered the inverse probability censoring weighted response ∆ i N i (t)/G(T i ∧ t|X i , Z i ) in the estimating equations and proposed to estimate the regression coefficients α(t) and γ by solving the estimating equations simultaneously. We denote the estimates asα(t) andγ. Under regularity conditions we showed that √ n(γ − γ) and √ n{α(t) − α(t)} are jointly asymptotically Gaussian and have the same limit distribution as respectively, where (G 1 , ..., G n ) are i.i.d. standard normals, τ is the of study time point, and explicit expressions forĈ γ ,Î α (t),Ŵ 1i (t) andŴ 2i (t) are given in .
For a given set values of covariates, (x, z), the predicted cumulative incidence function can be estimated byP 1 (t; x, z) = h −1 x α(t) + g(z,γ, t) , and we showed in SZG that √ n{P 1 (t; x, z)− P 1 (t; x, z)} has the same limit as where (G 1 , ..., G n ) are standard normals andŴ 3i (t; x, z) is a residual that can be estimated based on the data (see Scheike and Zhang 2008 for details).
Resampling techniques can be applied to construct (1−α)100% confidence bands for α j (t), j = 1, . . . , p, and P 1 (t; x, z), and to compute the p value of testing is an estimated standard error ofα(t).

Goodness-of-fit test
The FG model is commonly used for analyzing competing risks data and this model assumes that all covariates have constant effects over time (Beyersmann et al. 2009). Recently, we developed a goodness-of-fit test  for testing whether or not this is a reasonable assumption. We consider an extended version of the FG model where some effects are proportional as in the FG model (γ) and some effects are allowed to change their effects on the cumulative incidence function over time (α(t)). Therefore testing for example H 0 : α j (t) = β j , for all t ∈ [0, τ ] will determine whether the effect of x j is constant. Further, plotting the estimatedα j (t) with its confidence band will give a good idea about whether or not the proportionality assumption is satisfied or violated.
To test H 0 there are many possibilities, a simple test that relies only onα j (t) is to look at for j = 1, . . . , p. We derived the asymptotic distribution of this test process and proposed to compute the p value of the test based on a Kolmogorov-Smirnov type test-statistic sup t∈[0,τ ] |T j (t,α j )| or by a Cramer von Mises type test-statistic τ 0 {T j (s,α j )} 2 ds. The Cramer von Mises test is an alternative to the Kolmogorov-Smirnov test. Anderson (1962) showed that in the case of the two sample test, the Cramer von Mises test is more powerful than the Kolmogorov-Smirnov test. We compute both tests in the comp.risk() function.
In addition, we can plot the observed test process (5) and simulated test processes under the null hypothesis to visually examine whether a specific covariate has a time-varying effect.
All large sample properties and resampling techniques used for the test statistics are given in SZG and .

Worked example: Follicular cell lymphoma study
We consider the follicular cell lymphoma data from Pintilie (2007) where additional details also can be found. The data set can be downloaded from http://www.uhnres.utoronto.ca/ labs/hill/datasets/Pintilie/datasets/follic.txt, and consists of 541 patients with early disease stage follicular cell lymphoma (I or II) and treated with radiation alone (chemo = 0) or a combination treatment of radiation and chemotherapy (chemo = 1). Disease relapse or no response and death in remission are the two competing risks. The patients ages (age: mean = 57 and sd = 14) and haemoglobin levels (hgb: mean = 138 and sd = 15) were also recorded. The median follow-up time was 5.5 years.
We first estimate the nonparametric cumulative incidence curve using the timereg package and the cmprsk for comparison. We specify the event time and the censoring variable in timereg's comp.risk() function as Surv(dftime,cause == 0). The regression model contains only an intercept term (+ 1). The cause variable gives the causes associated with the different events. causeS = 1 specifies that we consider type 1 events, and the censoring code is given by the cens.code variable. The times at which the estimates are computed/based can be given by the argument times = times1, the default is to use all cause "1" time points that are numerically stable.
The cumulative incidence curve estimations based on the cmprsk's cuminc() function and the timereg's comp.risk() function are both identical to the product-limit estimator in the case without covariates (Figure 1 a and b). Figure 1 (a) shows the cumulative incidence curves for the two causes estimated by the cmprsk package. In Figure 1 (b) we show that the comp.risk() function can also be used to construct 95% confidence intervals (dotted lines) and 95% confidence bands (broken lines) based on resampling which is not available in the cuminc() function. The R packages etm (Allignol et al. 2011) and mstate (de Wreede et al. 2010 can also be used to compute the cumulative incidence curve with 95% confidence intervals, but they do not provide confidence bands. R> library("timereg") R> library("cmprsk") R> out1 <-comp.risk(Surv(dftime, cause == 0)~+ 1, data = fol, + cause, causeS = 1, n.sim = 5000, cens.code = 0, model = "additive") R> pout1 <-predict(out1, X = 1) R> group <-rep(1, nrow(fol)) R> fit <-cuminc(fol$dftime, cause, group, cencode = 0) R> par(mfrow = c(1, 2)) R> plot(fit,main = "cmprsk", xlab = "Years (a)") R> plot(pout1,xlim = c(0, 30), xlab = "Years (b)", main = "timereg", + uniform = 3, se = 2) Both the subdistribution hazard approach and the direct binomial modelling approach are based on an inverse probability of censoring weighting technique. When applying such weights it is crucial that the censoring weights are estimated without bias, otherwise the estimates of the cumulative incidence curve may also be biased. In this example, we find that the censoring distribution depends significantly on the covariates hgb, stage and chemo and is well described by Cox's regression model. The fit of the Cox model was validated by cumulative residuals, see Martinussen and Scheike (2006) for further details. As a consequence using a simple Kaplan-Meier estimate for the censoring weights may lead to severely biased estimates. We therefore add the option cens.model = "cox" in the function call, this uses all the covariates present in the competing risks model in the Cox model for the censoring weights. More generally it has been established that regression modelling for the inverse probability censoring weights can be used to improve the efficiency .
We now use prop in the model option to fit the model We first fit a general proportional model allowing all covariates to have time-varying effects.
Only the covariates x in model (6) are defined in the function call below. The covariates z in model (6) are specified by a const operator.

OUTPUT:
Competing risks Model The tests of significance based on the nonparametric tests show that stage and age are clearly significant, chemo is borderline significant (p = 0.056) and hgb is not significant (p = 0.889) in the fully nonparametric model. Plot options of sim.ci and score can be used to plot estimated regression coefficients α j (t) with its 95% confidence bands and to plot the observed test process for constant effects and simulated test processes under the null, respectively.
R> plot(outf, sim.ci = 2) R> plot(outf, score = 1) Figure 2 shows the time-varying covariate effects (α(t) of model (6)). It is evident that these effects are not constant over time, effects are considerably pronounced in the early timeperiod. The 95% pointwise confidence intervals, as well as 95% confidence bands (sim.ci=2 in the plot call, 2 for broken lines). Figure 3 shows the related test-processes for deciding whether the time-varying effects are significantly time-varying or whether H 0 : α j (t) = β j can be accepted. The summary of these graphs are given in the output, and we see that stage, age and chemo are clearly time-varying, and thus not consistent with the Fine-Gray model. The p values related to these plots are given in the above output, and we see that the Kolmogorov-Smirnov (supremum) test leads to p values of 0.068, 0.007, and 0.000, for stage, age and chemo, respectively. Similarly, the Cramer von Mises test statistics based on the same score processes are 0.001, 0.001, and 0.090, respectively. These test statistics are described in detail in Section 2. Note that the two different summaries of the test processes by the Kolmogorov-Smirnov and Cramer von Mises tests statistics are consistent with the figures, and the overall conclusion is that none of the three variables have proportional Cox type effects. In reality the command plot(outf, Figure 3: Observed test process (black line) and simulated test processes under the null (gray lines). score = 1) that produces Figure 3 also leads to a similar plot for the baseline, but we have only plotted the covariate components of the models. To plot, for example, the second covariate (after the intercept), stage, of the model we give the command plot(outf, score = 1, specific.comps = 2). We see that hgb is well described by a constant and we therefore consider the model with hgb having a constant effect and the remaining covariates having time-varying effects.
This final model is fitted with the call R> outf1 <-comp.risk(Surv(dftime, cause == 0)~stage + age + chemo + + const(hgb), data = fol, cause, causeS = 1, n.sim = 5000, cens.code = 0, + model = "prop", cens.model= "cox") R> summary(outf1) The covariate hgb has a constant effect over time withβ = 0.00195. Note that hgb is nonsignificant (p = 0.627), as in the nonparametric model (p = 0.889) as well as in the FG model (p = 0.534) where all effects are constant over time (see below). The covariates stage, age and chemo all have significantly time-varying effects, and the estimates of the effects of stage, age and chemo are very similar to those of the fully non-parametric model shown in Figure 2.
To make a comparison of the predictions based on the FG model we also fit this model: R> outfg <-comp.risk(Surv(dftime, cause == 0)~const(stage) + const(age) + + const(chemo) + const(hgb), data=fol, cause, causeS = 1, + n.sim = 5000, cens.code = 0, model = "prop", cens.model = "cox") R> summary(outfg) We note that the effect of hgb is almost equivalent with that based on the more appropriate model (shown above). But the estimate could be severely biased due to lack of fit of the other covariates in the model, and could thus misrepresent important features of the data.

Competing risks Model
Finally, we compare the prediction for the FG model with that of the semiparametric model that gives a more detailed description of the effects. We consider predictions for two different patients defined by the newdata assignment below. Patient type I: disease stage I (stage = 0), 40 years old and without chemotherapy treatment (chemo = 0), and patient type II: disease stage II (stage = 1), 60 years old and the radiation plus chemotherapy combination treatment (chemo = 1).
We plot the predictions without pointwise confidence intervals (se = 0) and without confidence bands (uniform = 0). The predictions shows in Figure 4  incidence curves of relapse for a type I and a type II patient are plotted in solid and dotted lines, respectively.  Figure 5 (b) compares the predictions for a type II patient. The broken lines around the two predictions represent the confidence band based on the flexible model. Figure 5 is produced by the following code.

Discussion
The flexible competing risks regression model for the cumulative incidence curves are implemented in the comp.risk() function in the timereg package for R. These models are useful for a detailed analysis of how covariate effects predicts the cumulative incidence, and allows for a time-varying effect of the covariates. This is particularly useful for examining the fit of simpler models where covariate effects are assumed constant. The goodness-of-fit procedure leads to an asymptotically justified p value. Another nice feature is that the comp.risk() can deal with cluster structure.
The predict() function yields predictions with confidence intervals as well as confidence bands which are useful for the researchers.