A measure of explained risk in the proportional hazards model

1. INTRODUCTION

The proportional hazards model is used with survival data to determine important prognostic factors and to assess subject-specific relative risk. The model with time independent covariates is specified as

where t represents survival time, λ(t|x) is the hazard function conditional on a set of covariates denoted by the vector x, and β₀ is the vector of regression coefficients that determines the relationship between the covariates and the risk of death. The factor λ₀(t) represents the hazard when each component of the covariate vector x equals zero, and exp[β₀^Tx] is the relative risk for a subject with covariate profile x.

The ability to separate patients by risk is important for understanding therapeutic options. Explained risk, within the context of the proportional hazards model, gauges the ability to delineate patient risk by their covariate profile. If, however, after accounting for known risk factors, there remains significant unexplained heterogeneity, then the model determination of patient risk will be suboptimal. Thus, the assessment of explained risk is an important but underutilized component in the development of risk models in survival analysis. One explanation for the infrequent application of explained risk in the proportional hazards model is the lack of a consensus measure.

2. PREVIOUS APPROACHES TO EXPLAINED RISK

The coefficient of determination R² is the standard measure of explained risk in the normal linear model with uncensored data. Although the proportional hazards model can be specified through the semiparametric linear transformation family

(1)

where Λ₀(t) is the baseline cumulative hazard function and ν is a standard extreme value random variable, there are two barriers to attaining an R² type measure with this model and censored data. First, the transformed scale of the survival time Λ₀(t) is unknown, and second, some survival times are unobserved due to censoring.

Korn and Simon (1990, 1991) used loss and risk functions to develop a framework for the construction of explained risk measures. The loss function L(T,γ) measures the closeness of the survival time T to a parameter γ, and the risk function R_L[f] = min_γE_f[L(T,γ)] is the expected loss minimized over the parameter of interest. Through these concepts, explained risk is defined as

where the null density f₀ represents the survival time density when the covariates are ignored, and f_I is the minimum risk density. Thus, explained risk measures the reduction in risk due to the incorporation of covariates. Note that due to the expectation, the estimation of explained risk is based solely on model estimates.

If the survival distribution is unknown, the expectation may be replaced by the sample average. Korn and Simon (1991) call this statistic explained residual variation. In comparison to explained risk, explained residual variation is a function of survival time and the model estimates. It provides direct evidence of predictive accuracy (Henderson, 1995) and is more robust to model misspecification (Rosthoj and Keiding, 2004). However, applying explained residual variation to survival data is problematic since not all survival times are observable due to censoring, and the estimate is a function of the maximum follow-up time of the study, making it difficult to compare measures across studies.

To adapt explained residual variation measures to censored survival times, proposals have been developed that incorporate inverse probability censoring weights (Graf and others, 1999; O'Quigley and Xu, 2001). Measures that incorporate these censoring weights remain dependent on the maximum follow-up time, and if the censoring time is not independent of the covariates, the inverse probability censoring weight requires a separate model for the censoring distribution as a function of the covariates. These issues have limited the application of inverse probability weights to explained residual variation measures for censored data.

In this work, an explained risk measure is developed using an expected loss function derived under the proportional hazards specification indicated in (1). Prior to presenting the proposed approach, 3 previous measures of explained risk in the proportional hazards model are reported.

An estimated explained risk measure for the proportional hazards model was derived by Kent and O'Quigley (1988). Their measure was based on the entropy loss function and the Kullback–Leibler information gain. Entropy is a measure of uncertainty of the conditional survival time distribution

where using (1), the transformed survival time s = logΛ₀(t) represents the log baseline cumulative hazard function, and f(s|x;β) is the conditional density of an extreme value random variable with location parameter β^Tx and scale parameter 1. The term h(x) denotes the marginal covariate density.

The Kullback–Leibler information gain

is a comparison of entropy measures for the proportional hazards model with and without covariates. Information gain that contrasts entropy when β = β₀ to β = 0 is also referred to as explained randomness (Kent and O'Quigley, 1988). Entropy for the covariate model is approximately 1.5772 provided the extreme value specification is correct. The proportional hazards assumption, however, cannot simultaneously hold for model (1) and any model containing a subset of covariates. Due to this non-nesting property, entropy for the null model requires simultaneous estimation of the intercept, slope, and the scale parameters from the extreme value regression model (Heinzl, 2000). The estimates of these parameters are derived as implicit solutions to score equations. Due to its lack of simplicity, Kent and O'Quigley (1988) propose an approximation to the estimate of explained risk based on the product moment correlation coefficient with a standard normal error variance

In this statistic, $\hat{\mathbf{\beta }}$ is the maximum partial likelihood estimate, and ${\hat{\Sigma }}_x$ is the estimated variance–covariance matrix of the covariate vector. The estimated explained risk measure R_KO² is simple to compute using any proportional hazards software.

A second explained risk statistic is a censored data version of the likelihood ratio statistic,

where l is the partial likelihood, and d is the number of failures. Kent (1983) proposed this measure for parametric models in the uncensored data case. O'Quigley and others (2005) noted that this statistic was sensitive to the rate of censoring and proposed R_LR², replacing the number of subjects with the number of failures in the exponent of the likelihood ratio statistic. There is, however, no formal justification for this substitution. The measure R_LR² is also easy to compute using standard proportional hazards software.

A third measure of explained risk was developed by Xu and O'Quigley (1999). To provide this estimate, additional notation is needed. Let C represent the underlying censoring times, independent of the survival time T and covariate X, with survival function G(t). For each subject, the time Y_i = min(T_i,C_i) is observed along with the censoring indicator δ_i = I[T_i ≤ C_i]. The at risk indicator for subject i at time t is denoted by ψ_i(t) = I[Y_i ≥ t]. The Xu and O'Quigley estimate is defined as $1-\text{exp}(-2\hat{\Gamma }(\hat{\beta }\left)\right)$, where

This measure is a function of inverse probability censoring weights and may be sensitive to values in the right tail of the survival distribution. In addition, if the censoring times are dependent on some of the covariates, then an additional model is needed to estimate the inverse probability censoring weights.

For these statistics, either a modification or approximation to the metric is needed to patch the problems due to censoring or the non-nesting property of the proportional hazards model. In this work, a measure of explained risk for the proportional hazards model will be developed that incorporates the non-nesting property of the proportional hazards model and is unaffected by censoring times that are independent of survival times conditional on the covariates.

3. EXPLAINED RISK IN THE PROPORTIONAL HAZARDS MODEL

In this proposal, the explained risk is composed of entropy loss functions derived from the full model, null model, and degenerate model. Applied to the proportional hazards model, this requires special consideration since the proportional hazards assumption will hold for only a single model, which is taken here as the full model. To derive the null model, the Korn and Simon (1990, 1991) paradigm is employed. The density under the null model is the average of conditional densities from the full proportional hazards model

Since the full model is assumed correct, the null density is constructed to indicate the effect of ignoring the covariates in the computation of the risk function. The entropy loss function is developed under the extreme value distribution with location parameter θ and scale parameter 1

which is chosen for its connection to the proportional hazards model, as portrayed in (1).

The proposed measure of explained risk is a function of the minimum expected entropy under the full model

the null model

and the degenerate (minimum entropy) model

Collecting terms, the explained risk is

Under the assumption that f(s|x;β) is the conditional density of an extreme value random variable with location β^Tx and scale parameter equal to 1, i.e. the proportional hazards specification is satisfied, and the covariate vector is centered around zero, then H_x = 1.5772, H_I = 1, and

The resulting measure is

Substitution of the maximum partial likelihood estimate $\hat{\mathbf{\beta }}$ provides the estimate of R²(β). Details for the entropy calculations are provided in the supplementary material available at Biostatistics online.

The degenerate model is used to construct the lower bound for the entropy under the extreme value loss function. The entropy for the full model is constant and is equal to one plus Euler's constant, which is rounded to 1.5772. Its constancy is derived from the monotone transformation in (1) that sets the error to an extreme value random variable with mean 0 and scale parameter equal to 1. The null model entropy is the explained risk lost due to ignoring the covariates.

Four properties of R²(β) are highlighted below.

1. The driver of this statistic is the subject-specific relative risk. Thus, an interpretation within the proportional hazards framework is a measure of explained relative risk.

2. When censoring is independent of survival conditional on the covariates, the estimate $R^2\left(\hat{\mathbf{\beta }}\right)$ is unaffected by the underlying censoring distribution and the length of the study. The estimate is obtained through the maximum partial likelihood estimator $\hat{\mathbf{\beta }}$ and the empirical distribution for the covariate vector. Since conditionally independent censoring has little effect on $\hat{\mathbf{\beta }}$, the censoring distribution has virtually no effect on the estimated explained risk measure.

3. R²(β) is bounded between 0 and 1. This result follows from an application of Young's inequality (Arnold, 1989)

where f is the convex function f(a) = exp(a) − 1 and g is its Legendre transform g(b) = 1 − b + blogb. Substituting a = β^Tx and b = 1, it follows that for each subject,

and since the {x_i} are centered around zero,

4. A device that is often cited to assess the adequacy of a newly developed R² measure is its comparability to the coefficient of determination in the normal linear model. When R² is applied to the normal linear scale transformation model

where m is an unknown monotone function and {ϵ_i} are independent, identically distributed, normal random variables with mean 0 and variance σ², then , and

This produces the measure

which for large n approximates the population value of the coefficient of determination for the normal linear scale transformation model

Thus, the proposed statistic is based on a primary component of the proportional hazards model, the relative risk function. It ranges between 0 and 1, permitting a standard metric for interpretability, and since it is stable in the presence of conditionally independent censoring, allows for comparability between studies. In addition, the simplicity of the estimate $R^2\left(\hat{\mathbf{\beta }}\right)$ enables a straightforward derivation of its asymptotic normal distribution, and ultimately, an asymptotic confidence interval for its population value ρ²(β) = lim_n→∞E[R²(β)].

THEOREM Assume the {x_i} are generated independent and identically distributed with E(x₁^Tx₁) < ∞ and μ_x(β) = E[exp(β^Tx₁)] > 0. Then

The estimate of the asymptotic variance of $R^2\left(\hat{\mathbf{\beta }}\right)$ follows by substituting the maximum partial likelihood estimate $\hat{\mathbf{\beta }}$ for β₀ in each component. The derivation of this asymptotic distribution is provided in the supplementary material available at Biostatistics online.

4. CONTRIBUTION OF A SUBSET OF FACTORS TO THE EXPLAINED RISK

In addition to the overall assessment of explained risk, the relative importance of a subset of factors from the full model may be ascertained. The full proportional hazards model is denoted by

and the density of a subset model, obtained by ignoring the heterogeneity due to the covariate vector z, is denoted by

The extreme value entropy loss function is now written as

where θ_x = β^Tx and θ_z = γ^Tz. The entropy for the full model is H_x,z = 1.5772, and the entropy for the subset model, ignoring the heterogeneity due to z, is

Therefore, the contribution of the covariate vector z to the explained relative risk

is equal to

Its estimate is obtained through substitution of the partial likelihood estimates $(\hat{\mathbf{\beta }},\hat{\mathbf{\gamma }})$ from the full model.

5. SIMULATIONS

A simulation study was performed to assess the properties of the explained relative risk measure and compare the results to the likelihood ratio, Kent and O'Quigley, and Xu and O'Quigley measures of explained risk described in Section 2.

The survival times were generated from the regression model

where the scalar covariate X_i was a standard normal random variable, and the ϵ_i were independent identically distributed Weibull random variables with scale parameter 1 and shape parameter taking on the values {1.20,0.80,0.50,0.30}. These parameters were chosen to create a wide range for the explained risk measures. The regression model produces a proportional hazards relationship between the covariate and survival time. Censoring times were generated independently of the covariate or survival time, from a uniform distribution with support (0,τ), where τ was varied to produce average censoring proportions approximately equal to {0,0.25,0.50,0.75}. For each simulation, the sample size was 100 and the average of 5000 replications was reported.

To determine the accuracy of these measures, the survival times were log transformed to produce a log-linear model, and the coefficient of determination

was computed in the uncensored data case, where $\stackrel{\sim }{\mathbf{\beta }}$ represents the least squares estimate. Under the Korn and Simon (1991) classification, R_LLM² is a measure of explained residual variation and is expected to be smaller, more pessimistic, than the explained risk statistics under study (Henderson, 1995). This statistic, which is common in the uncensored data case, does not account for censoring and hence the need for alternative measures.

The results in Table 1 indicate that the proposed explained relative risk R², and the Kent and O'Quigley measure R_KO² are unaffected by the percent censored, whereas the likelihood ratio R_LR², and Xu and O'Quigley measure R_XO² fluctuate as the censoring proportion (maximum follow-up time) varies. In the uncensored data case, the explained relative risk, the likelihood ratio statistic, and the Xu and O'Quigley measures were comparable to R_LLM² over the range of Weibull shape parameters explored. The likelihood ratio statistic, and the Xu and O'Quigley measures, however, increased as the percent censoring increases. The Kent and O'Quigley statistic R_KO² had a greater bias than the explained relative risk R², when compared to R_LLM² and R_LR² in the uncensored case. Since it is unaffected by censoring, the bias in R_KO² remained across censoring rates. These results indicate that the explained relative risk R² statistic is both stable and accurate relative to the other commonly cited measures of explained risk.

The estimated standard error of the proposed R² is comparable to the simulation standard error of the estimate. Its estimated variability (ASE) increased as the censoring proportion increased and as the Weibull shape parameter decreased. The latter result stems from the property that under the generated Weibull regression model, a decrease in the shape parameter indicates an increase in the variability of the underlying survival times.

Additional simulations (not shown) were generated for survival times from a proportional hazards model with a baseline cumulative hazard function equal to Λ₀(T) = log(1 + T). Thus, the survival times were not generated from a Weibull distribution, but the data still maintain the proportional hazards structure. The results were similar to the Weibull simulations in Table 1.

6. PROSTATE CANCER DATA EXAMPLE

A model was developed on 1006 castrate resistant metastatic prostate cancer patients, using 10 prognostic factors within a proportional hazards model (Armstrong and others, 2007). The model produced a concordance index (Harrell and others, 1984), a measure of model discrimination for the survival time, equal to 0.38 when transformed to a [0,1] scale. This result indicates that there remains important unknown or unrecorded factors that explain the heterogeneity in relative risk. In 2009, a study was undertaken to explore the role of a new biomarker, circulating tumor cells (CTC), as a prognostic factor for survival time in this prostate cancer population (Scher and others, 2009). CTC is a blood-based assay that provides information on the accumulation of tumor cells in the peripheral blood.

In the 2009 study, a proportional hazards model was generated from a cohort of 136 patients with castrate resistant metastatic prostate cancer. In addition to CTC, the factors considered for this model were: hemoglobin, prostate-specific antigen (PSA), lactate dehydrogenase (LDH), alkaline phosphatase, and albumin. A logarithmic transformation of some of the factors was used to account for their positively skewed distribution. The results from this analysis are presented in Table 2a. The test of the proportional hazards assumption developed by Grambsch and Therneau (1994) provided no evidence that the proportional hazards assumptions were violated for the full model (p = 0.681). The factors LDH (p < 0.001) and CTC (p = 0.004) were the strongest prognostic factors, indicating that CTC adds to understanding patient risk in the metastatic prostate cancer population. However, the magnitude of the p-value is sensitive to sample size, and it alone is not sufficient to evaluate a model's predictive power of risk assessment. To assess the model's adequacy in explaining relative risk, the evaluation provided in Table 2a is complemented with the calculation of the proposed R² statistic in Table 2b.

For the 6 factors under study, the explained relative risk is R² = 0.581(se = 0.086). The R-square value for the other approaches examined in the simulations were: R_LR² = 0.579, R_KO² = 0.540, and R_XO² = 0.607. Table 2b indicates the influence of the covariate subsets in the determination of R². As expected, the 2 factors that individually provide the greatest contribution to explained relative risk are LDH and CTC. The combination of LDH and CTC generated an R² = 0.411, indicating that over two-thirds of the model's explained relative risk was generated by these 2 factors. Adding hemoglobin to these 2 factors produced an R² = 0.547, representing 94% of the explained relative risk. The incremental value of the other factors (PSA, alkaline phosphatase, and albumin) to the explained relative risk was small.

7. DISCUSSION

The measure of explained relative risk for the Cox model is appealing due to its computational simplicity, and it is unaffected by censoring that is independent of survival conditional on the covariates. This stability stems from the proportional hazards specification and the fact that the relative risk is independent of follow-up time. The robustness of this measure may also be seen as a drawback, since the assumption of a constant relative risk with respect to time may be deemed too restrictive. From the analyst's standpoint, a prerequisite to its application is that goodness of fit methods such as Grambsch and Therneau (1994) or Lin and others (1993) have been applied and no violations of the model assumptions are found.

Royston (2006) states that the properties of a good explained variation measure are: (i) approximately independent of the amount of (independent) censoring; (ii) reducing to R² for the normal linear model; (iii) maintaining the nesting property for a subset of variable, i.e. model₁⊂model₂ then R₁² < R₂²; (iv) R² increases with the strength of the association; (v) construction of a confidence interval for the measure. With the exception of property 3, due to the non-nesting of proportional hazards models, the proposed R² measure satisfies these conditions.

Historically, patients with a given disease were classified into a small number of clinically defined stages that defined patient risk. In recent years, a greater understanding of the disease process has led to an increase in the granularity of risk assessment. It is, however, unclear whether this increase in model complexity has resulted in a greater accuracy in the determination of patient risk. If these models are accepted uncritically, it is problematic when they are used in clinical decision making, for example as selection criteria for entry into a clinical trial or as trial stratification variables. A highly accurate risk model is necessary for its confident application. The explained relative risk statistic, which measures the reduction in randomness due to the introduction of patient risk factors, provides a measure of accuracy in the presence of conditionally independent right censored data and is straightforward to apply.

SUPPLEMENTARY MATERIAL

Supplementary material is available at http://www.biostatistics.oxfordjournals.org

FUNDING

National Institute of Health (P50 CA92629).

Supplementary Material