A pseudo-R²measure for selecting genomic markers with crossing hazards functions

Notations

Let the random variables X and C be the failure and censoring times, and T = min(X, C) be the observed follow-up time. The random variables T and C are assumed to satisfy the condition of independent censoring [8]. We denote {N _i (t), t ≥ 0} the counting process that indicates the number of events that have occurred in the interval (0, t] for subject i, i = 1,⋯, n, so that N _i (t) takes values 0 or 1. Let Y _i be the at-risk process, so that Y _i (t) = 1 indicates that subject i is at risk just before time t, and Y _i (t) = 0 otherwise.

Let dN _i (t) = N _i (t ^- + dt) - N _i (t ^-) be the number of events occurring in the interval [t, t + dt) for subject i, the total number of events that have occurred in the interval (0, t] and the number of subjects at risk at time t. Finally, let represent the value of the g^th covariate for individual i.

Motivational situation: the modulating effect

In the following, we show how a simple interplay between two binary markers Z ⁽¹⁾ and Z ⁽²⁾ can lead to marginal crossing hazard functions.

The joint distribution of Z ⁽¹⁾ and Z ⁽²⁾ is defined by:

It is assumed that the hazard function of subject i with and is given by

(1)

where λ ₀(t) is an arbitrary unspecified baseline hazard function, and α and γ are unknown regression coefficients.

Assuming that model (1) is the true one, the consequences of omitting Z ⁽²⁾ on the formulation of the observed hazards ratio relative to Z ⁽¹⁾ is described below. Expressing model (1) in terms of the conditional survival function given leads to:

where S ₀(t) is the survival function corresponding to the baseline hazard function λ ₀(t). The survival function given ( ) follows directly from Bayes' theorem, and the hazard function given ( ) can be easily deduced as:

It is worth noting that this latter expression can be obtained as the expectation of (1) taken over Z ⁽²⁾ given the at risk process. Finally, the hazards ratio relative to the values and is given by:

(2)

It appears from this expression that hazards may cross over time. More precisely, it is shown in Additional File 1 that when α and γ are positive and assuming a balanced joint distribution for (Z ⁽¹⁾, Z ⁽²⁾), the hazards ratio inverts at a given time in (0; +∞). Obviously, such a time-dependence cannot be properly handled by using the proportional hazards model to analyze the data.

Semi-parametric model

The proposed model defines the survival function of subject i with covariate Z _i as follows

(3)

where λ ₀(t) is an unspecified baseline hazard function and β an unknown regression parameter. It is a particular case of a model that was proposed for handling hazards ratio that invert over time [9], and it corresponds to a semi-parametric generalization of the Weibull distribution. For subject i; i = 1, ⋯, n, the model (3) can be written in terms of the hazard function

(4)

where is the cumulative baseline cumulative hazard function.

In the simple case of a covariate Z taking values 0 or 1, the hazards ratio (HR) of the two groups corresponding to Z = 1 and Z = 0 is equal to . If β > 0, this function is increasing from 0 to +∞ and takes value 1 for . In this case, the risk of occurrence of the event is smaller in group 1 than in group 0 for 0 ≤ t < τ, and becomes greater when t > τ. If β < 0, the hazards ratio is decreasing from +∞ to 0 and takes value 1 for t = τ as calculated above. The risk of occurrence of the event is thus greater in group 1 than in group 0 when 0 < t ≤ τ and becomes smaller for t > τ.

Thus, as expected, model (4) allows hazards to cross over time. Note that the survival functions cross at a time larger than the crossing time of the hazards, and may not cross at a finite time.

At time t and for an individual i, i = 1, ⋯, n, the first derivative of the partial log-likelihood with respect to β is the score function:

The score function evaluated for β = 0 can be written at time t as

(5)

with ω(s) = 1 + log{Λ₀(s)}.

Pseudo-R²measure

The goal of this section is to propose a pseudo-R² index that can be interpreted in terms of percentage of separability between patients according to their survival times and marker measurements under the crossing hazards model (4). The approach used below is based on the score function (5). It extends the particular case that we considered in a former work [1] where we assumed the classical PH model. The main idea is to note that the score can be rewritten as a separability quantity between patients experiencing or not the event of interest. More precisely, the quantities U _i in (5) can be rewritten as

With .

From this expression, we show that, for a given covariate Z at time t, the U _i can be expressed as the weighted difference between the value of the covariate of the patient observed to experience the event of interest and the mean of the covariates of the group of patients observed to not experience the event.

An estimation of the U _i is given by

where , Λ₀(t _i ) is estimated by the left-continuous version of the Nelson's estimator [10, 11], and δ _i , i = 1, ⋯, n is the indicator of failure at time t _i .

For distributional reason, instead of the U _i , we introduce the so called robust scores W _i [2] which expressions are

(6)

Where

The W _i can be estimated by

(7)

The sum over i of the robust W _i , i = 1, ⋯, n is identical to the sum of the U _i . However, as for the Cox model, the W _i are independent, while the U _i are not.

Finally, the index is equal to the robust score statistic divided by the number of distinct uncensored failure times k:

The index D ₀ is interpreted in terms of percentage of separability over time between the event/non-event groups. Its calculation is easy as it does not require the estimation of the parameter β of the crossing hazards model. We can easily demonstrate that 0 ≤ D ₀ ≤ 1.

It is worth noting that the index D ₀ can be interpreted as a pseudo-R² measure. In the linear regression model, the R² (coefficient of determination) can be directly linked to likelihood-related quantities such as the Wald test, the likelihood ratio and the score statistics (see [12]). These formal relationships provide different ways to interpret the R². In the framework of non-linear models, statisticians have searched for a corresponding index and different pseudo-R² statistics have been proposed for censored data. Our proposed index is an extension of the definition of the R² for survival model with crossing hazards which relies on the score statistic.

Simulation Scheme

A simulation study was performed to describe the behavior of the proposed index, , in finite samples generated under the crossing hazards model (3), and to compare it to the behavior of other existing indices. Different situations were considered, corresponding to different covariate distributions, regression parameter values and sample sizes. The influence of various censoring distributions was also investigated. The indices that were compared to include our previous index derived from a PH model (i.e. calculated according to the same approach than with ω(t _i ) = 1, [1]) and other most usual indices: Allison's index [3], its modified version [5], Nagelkerke [4] and Xu and O'Quigley's [6] indices. All of them are designed for a PH model and are denoted , and , respectively. The different elements defining a configuration were the following. For a given subject, the distribution of the covariate Z included in model (3) was either discrete (Bernoulli ℬ(0.5)) or continuous (log-normal with mean 0 and variance 1/4, or uniform . These three distributions of Z were standardized to have the same variance. Two distributions for the survival time X were considered. The first one was defined by model (3) with Λ₀(t) = t. It is equivalent to a Weibull parametric model with scale parameter η = 1 and shape parameter α = exp(βZ). The second one correspond to a log-normal distribution with mean equal to 0 and standard deviation equal to e ^-βZ . These two distributions allowed to simulate the crossing hazards phenomenon.

The coefficient β was given a value such that e ^β = 1, or 2, or 3, or 4. It can be noticed that, in the case of a Bernoulli variable Z, the hazard functions corresponding to e ^β = 2, or 3, or 4 cross when the survival function of group Z = 1 equal 0.78, 0.82, 0.85 respectively. The censoring variable C _i was assumed to be independent from X _i given Z _i and distributed according to either a uniform or exponential C _i ~ℰ(γ) distribution. The parameters r and γ were calculated in order to yield an expected overall percentage of censoring p _c equal to 0%, 25% and 50%. The sample size n was taken equal to 50, 100 or 500. Data were generated as follows. For each subject i, i = 1, ⋯, n, a value of the covariate Z _i was generated. Given that value, a survival time X _i was generated according to a either Weibull distribution , or a log-normal distribution . The censoring variable C _i was independently generated and the observed follow-up time T _i was calculated as min(X _i , C _i ). For each configuration 1,000 independent replications were generated.

Simulation Results

Figures 1, 2 and 3 display the results of the simulations obtained for , , and for n = 100, for a Bernoulli, a uniform and a log-normal distribution, respectively, according to the values of e ^β and the percentage of censoring p _c . Figures 4, 5 and 6 give the results for n = 50. The results for n = 500 are given in Additional Files 2, 3 and 4.

As seen from Figures 1, 2, 3, 4, 5 and 6 and Additional Files 2, 3 and 4, when β = 0, i.e. in the absence of covariates, the different indices are close to zero for n = 50, 100 and 500. Among the six indices, only shows a mean value increasing regularly with β. The means of the five other indices do not appear to increase with β and remain below 0.05 even for the highest value of e ^β = 4. When β ≠ 0, the mean values of for the different sample sizes are fairly stable.

The standard errors of the six indices are small when β = 0. The standard errors of are larger than those of the other indices when β ≠ 0 and slightly decrease as β increases. As expected, the standard errors of the different indices decrease when n increases.

The mean value of does not appear to be sensitive to the censoring rate. However, the standard error of this index moderately increases as the percentage of censoring increases from 0% to 50%. In addition, the results obtained with a log-normal distribution of the survival time X on one hand (see Additional Files 5, 6 and 7 for the case n = 100), and with an exponential distribution of the censoring variable, on the other hand (not shown) are very similar, concerning and the other indices as well.

Application of the index on real data

In this section, we illustrate the use of the proposed index by selecting transcriptomic prognostic factors having a crossing effect in a lung cancer study. We compare the selection to the one obtained when relying on the index calculated under a proportional hazards model.

Dataset

This series is composed of 74 patients who underwent surgery at the Hôtel-Dieu Hospital (AP-HP, France) between August 2000 and February 2004 for stage IB (pT2N0) primary adenocarcinoma or large cell lung carcinoma of peripheral location [13]. Relapse-free survival was defined as the time from surgery until disease-related death, disease recurrence (either local or distant), or last follow-up examination. The median relapse-free survival time was 63.8 months. The two years relapse-free survival was 80.3%[71.2%, 90.5%], and the five years relapse-free survival was 59.3%[47.2%, 74.5%]. For each patient, we considered the gene expression measurments of 51,852 transcripts (obtained using Aymetrix HU133 Plus 2.0; Aymetrix, Santa Clara, CA, USA) located on the autosomal chromosomes.

Selection of the variables

The genes were ranked according to the value of either or . We decided to focus our attention on the first 200 top-ranked transcripts in both cases. The lowest separability for both indices were close to each other (29.8% for and 29.1% for ). Only a small proportion of transcripts (5%) was common to both lists.

We then examined the biological processes that were significantly over-represented in the two sublists using the PANTHER (Protein ANalysis THrough Evolutionary Relationships) classification system [14]. Results showed that the two indices allow selecting genes involved in different biological processes (See Additional File 8). For the cell-cycle process, selected 25 transcripts (significantly higher number than the 5% expected by chance), whereas selected only 16 transcripts (not significantly higher number than the 5% expected by chance). The two lists of genes involved in cell cycle are given in Additional File 9.

Among the 25 cell-cycle related transcripts selected according to , we discussed the behavior of two genes, namely FGFR2, MCL1, known to be involved in complex biological pathways. In particular, we examined whether the crossing phenomenon (observed om Figures 7.a and 7.b) could potentially be related to some effect modification of other genes.

The gene FGFR2 (fibroblast growth factor receptor 2) is known to be involved in various cancer types [15] and low gene expression measurements have been reported as linked to a shorter survival in lung cancer [16]. The analysis of FGFR2 gene expression taking into account FGF4 gene expression, which is one of its ligand, suggested a potential modulating effect between the two genes. In the following, we reported the hazards ratio (HR) computed under the Cox PH model for the four groups resulting from dichotomizing the two variables at the median. We also displayed the Kaplan-Meier curves on Figure 7.c. As seen on this latter, patients with low expression (below the median) of both FGFR2 and FGF4 have the worst prognosis (reference group). When FGFR2 is highly expressed (above the median) and FGF4 lowly expressed, the survival is not significantly improved (HR = 0.532 [0.202, 1.399]). However, the over-expression of FGF4 significantly improve the survival (HR = 0.329 [0.112, 0.967]). Finally, patients with a high expression of FGFR2 and FGF4 have the best prognosis (HR = 0.103 [0.021, 0.516]).

In the same way, we discussed the interaction between MCL1 and BCL2, two anti-apoptotic genes belonging to the BCL-2 gene family. Considered initially as oncogenes, the prognostic impact of BCL-2/MCL-1 for various type of cancer is debated due to their dual function on cell death and cell proliferation (for a review, [17]). The anti-apoptotic effect is associated with resistance to chemotherapy, leading an adverse prognostic role in some cancers such as leukemia or advanced ovarian tumors. In contrast, the anti-proliferative activity of MCL-1 and BCL-2 is associated with a favorable prognosis effect in some early carcinomas, such as lung adenocarcinoma [18]. Moreover, the combined analysis of MCL1 and BCL2 gene expressions indicated a potential modulating effect between them. As seen on Figure 7.d., in our lung cancer study of early lung adenocarcinomas treated by surgery alone, patients with low expression (below the median) of both genes MCL1 and BCL2 have the worst prognosis (reference group). When MCL1 is highly expressed (above the median) and BCL2 lowly expressed, the prognosis is not significantly improved (HR = 0.533[0.202, 1.4103]). On the contrary, patients with low expression of MCL1 and high expression of BCL2 have a significantly improved survival (HR = 0.296[0.091, 0.962]). Finally, the over-expression of both MCL1 and BCL2 gives the best prognosis (HR = 0.189[0.051, 0.700]).

In these two examples, we could hypothesize that the crossing effect observed in the marginal analysis of FGFR2 or MCL1 is related to some potential effect modification linked to FGF4 or BLC2, respectively. This hypothesis is consistent with the known biological activity of those genes. FGFR2 encodes for a receptor, which needs one of its ligand (i.e. FGF4 ) for activation and biological activity. Also, BCL2 and MCL-1 encode two proteins of the same family, which may act together on the cell through heterodimerization on apoptosis or cell proliferation. The resulting subgroup of patients defined by high expression of these genes couples might be clinically relevant and the object of further investigations.