In this section, we first present a simple situation which motivates the use of the semi-parametric non-proportional hazards model introduced in the next subsection.
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 gth covariate for individual i.
Motivational situation: the modulating effect
In the following, we show how a simple interplay between two binary markers Z
(1) and Z
(2) can lead to marginal crossing hazard functions.
The joint distribution of Z
(1) and Z
(2) is defined by:
It is assumed that the hazard function of subject i with
and
is given by
where λ
0(t) is an arbitrary unspecified baseline hazard function, and α and γ are unknown regression coefficients.
Model (1) describes a modulating effect of the two markers Z
(1) and Z
(2), whereby Z
(2) has a multiplicative effect on the hazard and Z
(1) has a multiplicative effect only if Z
(2) equals one (so called effect modification). The corresponding hazard functions according to the values of Z
(1) and Z
(2) are shown in Table 1.
Assuming that model (1) is the true one, the consequences of omitting Z
(2) on the formulation of the observed hazards ratio relative to Z
(1) is described below. Expressing model (1) in terms of the conditional survival function given
leads to:
where S
0(t) is the survival function corresponding to the baseline hazard function λ
0(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
(2) given the at risk process. Finally, the hazards ratio relative to the values
and
is given by:
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
(1), Z
(2)), 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
where λ
0(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
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
with ω(s) = 1 + log{Λ0(s)}.
Pseudo-R2measure
The goal of this section is to propose a pseudo-R2 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
, Λ0(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
Where
The W
i
can be estimated by
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
0 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
0 ≤ 1.
It is worth noting that the index D
0 can be interpreted as a pseudo-R2 measure. In the linear regression model, the R2 (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 R2. In the framework of non-linear models, statisticians have searched for a corresponding index and different pseudo-R2 statistics have been proposed for censored data. Our proposed index is an extension of the definition of the R2 for survival model with crossing hazards which relies on the score statistic.