Evaluating prognostic biomarkers for survival outcomes subject to informative censoring

1. Introduction

Biomarkers play increasingly important roles in clinical research and practice. A biomarker is said to be prognostic if it is associated with clinical outcomes in a defined population under specified conditions. When the outcome of interest is a failure time, prognostic biomarkers are commonly evaluated using various concordance measures^1–4 and time-dependent receiver operating characteristic (ROC) curves.^3,5 A succinct summary of an ROC curve is the area under the curve (AUC).^6–8 In this article, we consider two specific performance measures: a C-index similar to that defined by Heagerty and Zheng ³ with an ability to handle discreteness in the distribution of marker values, and a time-dependent AUC where cases (controls) are defined as subjects who have (not) experienced a failure by a specified time point. In the terminology of Heagerty and Zheng, ³ the AUC we consider is a cumulative/dynamic AUC in the sense that the set of cases is cumulative and the set of controls is dynamic. We choose to focus on these performance measures because they are easy to interpret and commonly used in practice. In the rest of this article, the terms C-index and AUC refer specifically to these performance measures unless otherwise stated.

A major challenge in estimating these performance measures is that failure times are frequently subject to right censoring. Several methods have been developed to adjust for censoring in estimating the C-index^3,4 and the AUC.^9–12 These methods generally assume that censoring is uninformative in the sense that the censoring time is independent of the failure time with or without conditioning on the value of the prognostic marker being evaluated. The assumption of uninformative censoring is convenient to use but often questionable in clinical studies, where subjects may be lost to follow-up for various reasons. In such studies, appropriate handling of censoring frequently requires adjusting for additional covariates measured at baseline (i.e. start of follow-up).

As a motivating example, consider the expression level of the estrogen receptor gene estrogen receptor 1 (ESR1) as a prognostic biomarker for survival outcomes in breast cancer patients. ¹³ The prognostic value of ESR1 expression level in breast cancer has been studied extensively in terms of hazard ratios.^14–17 It is of interest to estimate the C-index and the AUC for this biomarker, and we aim to provide such an evaluation based on a cohort of 272 young ( < 53 years of age) women, with primary breast cancer of stage I or II, who were diagnosed and treated at the Netherlands Cancer Institute. ¹⁸ All patients were treated by modified radical mastectomy or breast-conserving surgery, including dissection of the axillary lymph nodes, followed by radiotherapy if indicated. ESR1 expression level was measured retrospectively on fresh-frozen tumor tissues. Follow-up information was extracted from the medical registry of the Netherlands Cancer Institute, with a median follow-up of 6.7 years. To assess the plausibility of the informative censoring assumption, the censoring time was regressed on important baseline covariates (see Table 1 of van de Vijver et al. ¹⁸ ) in addition to ESR1 under a Cox regression model. ¹⁹ A significant association was found between the censoring time and tumor diameter, suggesting that censoring was likely informative of disease characteristics and patient prognosis.

In this article, we develop and compare new methods for estimating the C-index and the AUC under informative censoring. We assume that the censoring time is conditionally independent of the failure time after adjusting for a collection of baseline covariates including the biomarker under investigation. Under this assumption, it is straightforward to obtain plug-in (PI) estimators based on a parametric or semiparametric model for the conditional distribution of the failure time given baseline covariates. This approach generalizes the AUC estimator of Chambless and Diao ⁹ based on a proportional hazards model from uninformative censoring to informative censoring. It is also straightforward to obtain inverse probability weighted (IPW) estimators based on a regression model for the censoring time conditional on baseline covariates. These IPW estimators are similar to those of Uno et al., ⁴ Hung and Chiang, ¹⁰ and Song et al. ¹² but are consistent under informative censoring. Motivated by semiparametric theory, we develop doubly robust (DR), locally efficient estimators that involve both a failure time regression model and a censoring time regression model, are consistent and asymptotically normal if either model is correct, and attain the nonparametric information bound if both models are correct. Though designed for informative censoring, the DR estimators remain applicable when the censoring is uninformative and provide robustness and efficiency improvements over existing methods for uninformative censoring.

The rest of the article is organized as follows. Section 2 defines the notations and states key assumptions. Sections 3 to 5 describe the PI, IPW, and DR methods, respectively. Section 6 reports a simulation study, and Section 7 presents the applications to real data. The article ends with a discussion in Section 8.

For a generic subject, let T be a failure time of interest (e.g. time to recurrence or death, whichever comes first) and X a quantitative biomarker (e.g. ESR1 expression level) intended to make predictions about T . Suppose that higher values of X are hypothesized to predict shorter times to failure; for ESR1, this condition can be met by negating the original marker values. In this notation, the C-index considered by Heagerty and Zheng ³ and Uno et al. ⁴ is P ( X 1 > X 2 | T 1 < T 2 , T 1 < τ ) , where the subscripts 1 and 2 denote two independent subjects and τ > 0 is a time point chosen to ensure identifiability. Since the distribution of X may have discrete components, we consider a slightly refined C-index given by

Estimation of θ c and θ a is complicated by the fact that T may be right-censored by a censoring time C . We assume that T and C are conditionally independent given W , written

A general motivation for PI estimators is provided by the following identity: for any real-valued function g such that E | g ( X 1 , T 1 , X 2 , T 2 ) | < ∞ , it follows from the law of iterated expectations that

In the PI approach, estimation of θ c and θ a requires estimating S T ( t | W ) , t ∈ ( 0 , τ ] . To deal with the curse of dimensionality, one may assume a parametric or semiparametric model for S T ( t | W ) , say S T ( t | W ; β ) , with a parameter β which may be finite- or infinite-dimensional. For example, S T ( t | W ; β ) may be specified as a parametric Weibull regression model or a semiparametric Cox regression model. Both models can be written as

Let β ^ be an estimate of β , which may be obtained from maximum likelihood or another estimation method. Then θ c can be estimated by

4. IPW estimators

It is straightforward to adapt existing IPW estimators of θ c and θ a ,^4,10,12 originally developed for uninformative censoring, to the present setting of informative censoring. Before describing specific estimators, we first give a general and intuitive motivation for such estimators. Note that both θ c and θ a are conditional expectations involving pairs of subjects and can be estimated using U-statistics in the absence of censoring. In the presence of censoring, we can start with an eligible pair of subjects with “full data” ( W 1 , T 1 , W 2 , T 2 ) and ask how their eligibility can be verified using the observed data ( W 1 , D 1 , U 1 , W 2 , D 2 , U 2 ) . Not all eligible pairs are verifiable, and only the verified eligible pairs enter into the IPW estimator. To adjust for a possible selection bias in the verification process, we weight each verified eligible pair by the inverse of their verification probability based on the full data.

In the definition of θ c , the conditioning event is T 1 < T 2 ∧ τ . This can be verified if and only if D 1 = 1 and U 1 < U 2 ∧ τ . In other words, the pair will be selected for IPW estimation if and only if

In reality, S C ( t | W ) is usually unknown but can be estimated using standard failure time regression models, treating C as the event time of interest and T as a censoring time for C . Suppose S C ( t | W ) is modeled using a parametric or semiparametric model, say S C ( t | W ; α ) , where α is a finite- or infinite-dimensional parameter. For example, S C ( t | W ; α ) may be specified as a Weibull regression model or a Cox regression model with

5. DR estimators

As noted in Section 3, each of the four quantities ( ν c , δ c , ν a , δ a ) that together determine ( θ c , θ a ) can be written as E { g ( X 1 , T 1 , X 2 , T 2 ) } for some function g . Specifically, we have g ( X 1 , T 1 , X 2 , T 2 ) = h ( X 1 , X 2 ) I ( T 1 < T 2 ∧ τ ) for ν c , I ( T 1 < T 2 ∧ τ ) for δ c , h ( X 1 , X 2 ) I ( T 1 ≤ τ < T 2 ) for ν a , and I ( T 1 ≤ τ < T 2 ) for δ a . In what follows, we consider DR estimation of ϑ = E { g ( X 1 , T 1 , X 2 , T 2 ) } for a generic function g before deriving DR estimators of ( θ c , θ a ) .

DR estimators are usually motivated by semiparametric theory.^20,21 In the present setting, the efficient influence function for estimating ϑ is given by

A general approach to DR estimation is to solve an estimating equation motivated by the efficient influence function. Specifically, we can estimate ϑ by solving the equation ∑ i = 1 n ϕ ^ eff ( O i ) = 0 , where ϕ ^ eff is an “estimate” of ϕ eff with unknown quantities (other than ϑ itself) replaced by suitable estimates. Suppose S T ( t | W ) is estimated by S T ( t | W ; β ^ ) as in Section 3, and S C ( t | W ) is estimated by S C ( t | W ; α ^ ) as in Section 4. Write Λ C ( t | W ; α ^ ) = − log { S C ( t | W ; α ^ ) } and M C ( t | W ; α ^ ) = ( 1 − D ) I ( U ≤ t ) − Λ C ( U ∧ t | W ; α ^ ) . There are different ways to estimate the function g ave , with different statistical properties. There is a PI estimator given by

Finally, let θ ^ c DR = ν ^ c DR / δ ^ c DR and θ ^ a DR = ν ^ a DR / δ ^ a DR , where the estimators ( ν ^ c DR , δ ^ c DR , ν ^ a DR , δ ^ a DR ) are obtained from equation (3) with appropriate choices for g . Using the continuous mapping theorem and the delta method, it is easy to see that the estimators ( θ ^ c DR , θ ^ a DR ) are consistent and asymptotically linear under the union model and nonparametrically efficient under the intersection model. Because the DR estimators involve multiple layers of integration, they can be quite computation-intensive, especially when the working models are such that integrals have to be evaluated numerically.

Here we report a simulation study of the finite-sample performance of the estimation methods described in Sections 3 to 5. Data are generated according to the following mechanism:

We take X = W 1 as the biomarker of interest and, in defining ( θ c , θ a ) , consider τ = t 0.2 and t 0.5 , where t p denotes the p th quantile of the marginal distribution of T . The corresponding true values of ( θ c , θ a ) range from 0.69 to 0.75. These parameters are estimated using the methods in Sections 3 to 5 with correct and incorrect working models for the conditional distributions of T and C given W . We consider both parametric Weibull regression models and semiparametric Cox regression models. The correct Weibull model for T is given by (4) with the same covariate vector W T and with λ T 0 ( t ) = κ T t ρ T − 1 , where κ T and ρ T are unknown parameters. The incorrect Weibull model for T has W T = ( W 1 , I ( W 2 > 0 ) ) ′ and is otherwise identical to the correct one. The correct Cox model for T is given by (4) with the same covariate vector W T and an unspecified baseline hazard function. The incorrect Cox model for T has W T = ( W 1 , I ( W 2 > 0 ) ) ′ and is otherwise identical to the correct one. The (correct and incorrect) models for C are specified in the exact same fashion as the models for T . All models are estimated using maximum likelihood.

The simulation results are summarized in Tables 1 (for θ c ) and 2 (for θ a ) in terms of empirical bias and standard deviation. When the working models are correct, the PI estimator is minimally biased, while the IPW estimator sometimes exhibits a large bias, which may be related to large variability. The DR estimator based on Cox models also shows visible bias at n = 200 , which appears to decrease with increasing sample size. When the working models are misspecified, the PI and IPW estimators become more biased, while the DR estimator is more protected when one of the working models is misspecified. Even when both models are misspecified, the DR estimator is frequently less biased than the PI and IPW estimators. Regardless of model correctness, the PI estimator has less variability than the IPW and DR estimators. Theoretically, the DR estimator is expected to be more efficient than the IPW estimator, at least when the working models are correct. This efficiency advantage is shown in Table 2 but not in Table 1, possibly due to limited precision of numerical integration in the DR estimator. For any method, variability generally increases when the censoring rate is higher and when the working models are specified as Cox regression models instead of Weibull regression models.

OPEN IN VIEWER

Table 2.

Simulation results for estimating the AUC ( θ a ): empirical bias and SD multiplied by 100.

			P ( T > C ) ≈ 0.57				P ( T > C ) ≈ 0.38
Estimation	Model	Model	τ = t 0.2		τ = t 0.5		τ = t 0.2		τ = t 0.5
method	for T	for C	Bias	SD	Bias	SD	Bias	SD	Bias	SD
Weibull regression, n = 200
PI	Correct		0.26	3.10	−0.01	3.53	0.21	2.70	0.00	2.88
PI	Incorrect		0.99	3.03	3.05	2.94	0.03	2.67	2.13	2.56
IPW		Correct	0.50	4.93	2.93	6.65	0.37	4.43	0.84	4.69
IPW		Incorrect	1.34	4.56	3.22	6.29	0.87	4.34	1.80	4.21
DR	Correct	Correct	0.12	4.61	−0.17	5.99	0.25	4.37	−0.05	4.44
DR	Correct	Incorrect	0.14	4.58	0.13	5.47	0.25	4.37	0.01	4.07
DR	Incorrect	Correct	0.40	4.68	1.35	5.97	0.37	4.38	0.53	4.47
DR	Incorrect	Incorrect	1.03	4.51	2.45	5.30	0.72	4.33	1.44	3.96
Weibull regression, n = 500
PI	Correct		0.29	1.94	0.20	2.24	0.23	1.68	0.15	1.81
PI	Incorrect		1.03	1.91	3.11	1.91	0.06	1.68	2.16	1.65
IPW		Correct	0.29	3.29	2.59	5.19	0.09	2.90	0.78	3.44
IPW		Incorrect	1.31	3.01	3.32	4.32	0.67	2.84	1.91	2.74
DR	Correct	Correct	−0.03	3.11	−0.19	4.03	0.06	2.89	−0.02	3.09
DR	Correct	Incorrect	0.02	3.07	0.28	3.55	0.07	2.87	0.12	2.60
DR	Incorrect	Correct	0.00	3.35	1.14	3.99	0.09	2.88	0.41	3.05
DR	Incorrect	Incorrect	0.86	3.00	2.52	3.50	0.49	2.84	1.50	2.55
Cox regression, n = 200
PI	Correct		0.19	3.18	−0.06	3.66	0.15	2.77	−0.06	2.98
PI	Incorrect		1.00	3.11	3.02	2.99	0.10	2.76	2.13	2.63
IPW		Correct	0.58	4.89	3.14	6.35	0.39	4.41	0.97	4.55
IPW		Incorrect	1.37	4.56	3.42	6.08	0.86	4.34	1.88	4.15
DR	Correct	Correct	0.92	4.74	1.28	5.78	0.58	4.36	0.58	4.24
DR	Correct	Incorrect	0.84	4.69	1.32	5.36	0.51	4.35	0.50	4.02
DR	Incorrect	Correct	1.22	4.78	2.60	5.49	0.71	4.36	1.17	4.19
DR	Incorrect	Incorrect	1.71	4.61	3.55	5.10	0.97	4.31	1.92	3.90
Cox regression, n = 500
PI	Correct		0.26	1.97	0.17	2.29	0.19	1.70	0.11	1.85
PI	Incorrect		1.08	1.94	3.13	1.95	0.15	1.71	2.21	1.67
IPW		Correct	0.30	3.29	2.71	4.94	0.10	2.89	0.81	3.41
IPW		Incorrect	1.31	3.01	3.42	4.22	0.66	2.83	1.93	2.69
DR	Correct	Correct	0.43	3.05	0.65	4.05	0.21	2.88	0.30	2.97
DR	Correct	Incorrect	0.49	3.16	1.04	3.79	0.19	2.89	0.34	2.60
DR	Incorrect	Correct	0.51	3.11	1.88	3.97	0.24	2.88	0.75	3.09
DR	Incorrect	Incorrect	1.31	3.08	3.25	3.65	0.60	2.86	1.72	2.56

AUC: area under receiver operating characteristic curve; SD: standard deviation; PI: plug-in; IPW: inverse probability weighting; DR: doubly robust.

We now illustrate the proposed methods with real data from two studies. The datasets and the R code for data analysis are available at https://github.com/zhi-wei-zhang/ProgSurv.

This article provides methods to account for informative censoring in prognostic biomarker studies with right-censored survival outcomes. The PI and IPW methods are straightforward extensions of existing methods for uninformative censoring. The DR method appears new to this literature (on prognostic biomarker evaluation). While designed for informative censoring, the DR method is obviously applicable to the setting of uninformative censoring. Based on theoretical considerations and simulation results, the IPW method does not have a clear advantage, and the PI and DR methods appear more promising. An important consideration in choosing between the PI and DR methods is the relative importance of bias versus variance. If one is primarily concerned about bias, the DR method is worth recommending because of its extra robustness. Even with the DR method, one cannot be too careful in specifying the working models, as the simulation results demonstrate that the DR method can be seriously biased when both working models are severely misspecified. On the other hand, the PI method does have an efficiency advantage when the failure time regression model is correctly specified, and may be preferable if one is more concerned about efficiency or mean squared error. The PI method is also simpler, easier to explain and implement, and computationally less expensive than the DR method. These practical issues, together with the bias-variance trade-off, should be considered carefully in choosing an estimation method.

Informative censoring may be a lesser issue in clinical trials with negligible drop-out, where the observation window is largely determined by the subject-specific enrollment date and the end date of study follow-up. Informative censoring is more likely to occur in observational studies where drop-out decisions may be influenced by individual health conditions. As illustrated in Section 1, a simple way to evaluate the extent of informative censoring is to regress the censoring time on a set of baseline covariates that are known or suspected to be prognostic for the medical event of interest. Significant associations found in this regression analysis would indicate the presence of informative censoring. The lack of significant associations does not guarantee the absence of informative censoring but may help to alleviate concerns about informative censoring.

If one decides to account for informative censoring induced by baseline covariates, one then has to choose a set of baseline covariates to adjust for. The baseline covariates that require adjustment are those that are plausibly associated with both the medical event of interest and the censoring time. In clinical studies, there is usually a set of prognostic variables with known or suspected associations with clinical outcomes; these are typically summarized in Table 1 of a medical research report. Some of these may be associated with censoring, though such associations may be less well understood because censoring itself is not of scientific interest. There may be insufficient data and background information to accurately identify a subset of prognostic variables that are also associated with censoring. In that case, we recommend that all prognostic variables be included in the adjustment for informative censoring if the number of such variables is not too large. We do not recommend the use of a data-driven variable selection procedure until its statistical properties are well understood, theoretically and empirically.

The main limitation of the proposed methods is their reliance on the conditional independence Assumption (1). This assumption cannot be validated with observed data and must be based on substantive knowledge. Further research is needed to develop sensitivity analysis methods that address possible violations of Assumption (1). The DR method is further limited by computational difficulties due to numerical integration. How to ease the computational burden of the DR method is another question for future research. In this article, methods are implemented using working models that converge at the parametric rate. In related areas (e.g. causal inference), there is growing interest in the use of nonparametric machine learning methods which provide greater flexibility and variable selection capability at the expense of slower convergence.^23–26 It will be of interest to adapt some of those approaches to the present setting of prognostic biomarker evaluation.

Supplemental Material