In practice, some censoring indicators may be unobserved, and we develop methods for this situation. For example, we analyzed data from a clinical trial involving patients with glioblastoma, a form of brain cancer. One measure of post-treatment quality of life is the length of time until a patient declines to a state in which mobility is lost and the cancer has returned. We evaluated data on adult glioblastoma patients over 50 years of age. Progression status was always known, but ambulatory status was unknown in some cases. By the end of the study, some patients were non-ambulatory and had progressed, some had progressed but remained ambulatory, some had progressed but their ambulatory status was unrecorded, and some had not progressed. Our analysis focused on the time to non-ambulatory progression, and these four types of events contributed uncensored observations, censored observations, observations with a missing censoring indicator, and censored observations, respectively. In particular, we were interested in how sex and age affected the time to non-ambulatory progression. Our methods have applications in many other areas, including analyses with missing death certificates in epidemiological studies, unperformed necropsies in animal experiments, and unidentified component failures in quality control settings.

The methods we propose are formulated in the context of a censored survival problem with missing censoring indicators, but they can be used to analyze competing risks data on event-specific failure times when some event indicators are unknown. In practice, the usual competing risks analysis focuses on a particular event type of interest and treats the times to failure from other events as censored observations on the event time of interest. Thus, one can view an unknown failure type as a missing censoring indicator and use our methods.

If censoring indicators are missing, standard estimators for censored data are not directly applicable. One solution is to ignore cases with missing data and apply conventional analyses to the complete cases only. However, the efficiency of this complete-case approach decreases as the degree of missingness increases. Also, the complete-case estimator is consistent only when censoring indicators are missing completely at random (MCAR). Alternatively, some analyses directly incorporate survival data with missing censoring indicators. For example, Dinse (1982); Lo (1991); McKeague and Subramanian (1998), and Subramanian (2004, 2006) proposed survival estimators. Within the regression context, Goetghebeur and Ryan (1990, 1995); Dewanji (1992); Lu and Tsiatis (2001), and Tsiatis et al. (2002) developed methods under a proportional hazards model; Zhou and Sun (2003) and Lu and Liang (2008) worked with an additive hazards model; and Gao and Tsiatis (2005) employed a linear transformation competing risks model. Linear regression techniques, however, have received little attention in situations where censoring indicators are missing, and thus we now focus on that particular problem.

Suppose there are n subjects and their observations follow the linear regression model where Y_i is a response variable, X_i is a vector of d covariates, β is a d-vector of regression coefficients, and ε_i is a random error for subject i(i = 1,…,n). Given {X₁,…,X_n}, the errors {ε₁,…,ε_n} are independent and identically distributed (i.i.d.) with mean zero. Let {C₁,…,C_n} be i.i.d. censoring variables, where C_i can randomly censor Y_i on the right. The observed response is Ỹ_i = min(Y_i, C_i), the censoring indicator is δ_i = I (Y_i ≤ C_i), and the missingness indicator ξ_i is 1 if δ_i is observed and 0 otherwise. Although calling δ_i a censoring indicator might seem backwards, as δ_i = 1 indicates an uncensored observation rather than a censored observation, we follow the convention that refers to δ_i as a censoring indicator.

We observe {Ỹ, X, δ} for subjects with complete data (ξ = 1) and {Ỹ, X} for subjects with a missing censoring indicator (ξ = 0). Typically the missingness mechanism is assumed to produce data that are missing at random (MAR), which implies pr (ξ = 1|Ỹ, X, δ) = pr (ξ = 1|Ỹ, X). Alternatively, we assume pr (ξ = 1|Ỹ, X, δ) = pr (ξ = 1|Ỹ), which is more restrictive than the MAR assumption but less restrictive than the common MCAR assumption that pr (ξ = 1|Ỹ, X, δ) = pr (ξ = 1). We made an assumption stricter than MAR to avoid postulating a parametric model for pr (ξ = 1|Ỹ, X) and to sidestep problems encountered when applying a nonparametric method to sparse data, such as can arise in high-dimensional situations. Under the weaker MAR assumption, we expect to obtain similar results for a reasonable model in the parametric case or with sufficient data in the nonparametric case, but we leave such analyses for future work.

In this paper, we propose three estimators for β. Our methods extend an approach taken by Koul et al. (1981) when all censoring indicators are observed. They used the formula to transform observed data (Ỹ_i, δ_i) into synthetic data , where G_n(c) is a product-limit type estimator of the censoring distribution G(c) = pr(C ≤ c). Koul et al. estimated β by applying a modified least squares procedure to , but their estimator does not handle missing censoring indicators. Also, they deflate the synthetic response to 0 when the true response Y_i exceeds C_i, and they inflate beyond Y_i when Y_i is actually observed, both of which are counterintuitive. Alternatively, our estimators ease this dilemma by replacing δ_i in (1) with a weighted average of δ_i and its estimated conditional mean _i, say p_iδ_i + (1 − p_i)_i, which allows statistical inference about β when some censoring indicators are missing.

For known m(·) and G(·), this result suggests the following least squares estimator of β:

In practice, however, m(·) and G(·) are usually unknown, and hence _R is not defined. One simple solution is to calculate _R using estimates of m(·) and G(·).

For example, consider a parametric model for m(z) of the form m₀(z, θ), where m₀(·, ·) is a known continuous function and θ is a vector of unknown parameters. The parameter vector θ can be estimated by _n, the maximizer of the likelihood and thus m(z) can be estimated by _n(z) = m₀(z, _n). As for estimating G(·), we begin by defining π(ỹ) = E (ξ|Ỹ = ỹ), which can be estimated nonparametrically by where W(·) is a kernel function and b_n is a bandwidth sequence. Next, we define a Horvitz and Thompson (1952) type inverse probability weighted estimator of u(ỹ) = E(δ|Ỹ = ỹ): where K (·) is a kernel function and h_n is a bandwidth sequence. Finally, similar to Dikta (1998) and Wang and Ng (2008), we define the following estimator of G(ỹ): where R_i denotes the rank of Ỹ_i (i = 1,…,n). Now, analogous to _R, we can define the following least squares regression calibration estimator of β: where Ŷ_i,R is Y_i,R with m(·) and G(·) replaced by _n(·) and Ĝ_n(·), respectively.

Define for 1 ≤ r, s ≤ d + 1 and assume:

The following theorem addresses the asymptotic normality of _R.

Theorem 2.1 Under assumptions (C.1), (C.2), (C.3) and (C.4), we have where V_R = Σ⁻¹Ω_RΣ⁻¹, Σ = E(XX), and Ω_R is defined in (A.15) of the Appendix.

The structure of Ω_R is complex; thus, estimating V_R by replacing the unknowns in Ω_R with appropriate component estimates is difficult to implement in practice. Alternatively, we estimate V_R by n_J, where _J is a jackknife estimator of the asymptotic variance of _R. In particular, we use a jackknife variance estimator of the form proposed by Peddada and Patwardhan (1992), which was motivated by a minimum norm quadratic unbiased estimator (MINQUE) and which in our case reduces to the simple form:

Theorem 2.2 Under assumptions (C.1), (C.2), (C.3) and (C.4), we have

Analogous to Sect. 2, we define a least squares imputation estimator for β:

Under the assumed missingness mechanism, _I also can be motivated by the relationship where Y_I = {ξδ + (1 − ξ)m(Z)}Ỹ/{1 − G(Ỹ)} and Y_R = m(Z)Ỹ/{1 − G(Ỹ)}.

The following theorem addresses the asymptotic normality of _I.

Theorem 3.1 Under assumptions (C.1), (C.2), (C.3) and (C.4), we have where V_I = Σ⁻¹Ω_IΣ⁻¹, , and is defined in (A.20) of the Appendix.

In practice, m(·), π(·) and G(·) are usually unknown, but we can use estimates of these quantities to calculate inverse probability weighted synthetic data:

Note that Ŷ_i,W is very similar to Ŷ_i,I, except the former weights each ξ_i inversely proportional to the estimated conditional mean of ξ_i (i.e., the estimated conditional probability that δ_i is known, given Ỹ_i). Applying least squares techniques to these synthetic data produces the inverse probability weighted estimator

The following theorem addresses the asymptotic normality of _W.

Theorem 4.1 Under assumptions (C.1) through (C.7), we have where V_W = Σ⁻¹Ω_WΣ⁻¹, , and is defined in (A.26) of the Appendix.

As described earlier, V_W can be estimated by jackknife methods. We see _W has a larger asymptotic variance than _I, and hence _R, by noting for any d-vector α. One explanation is that if Y_i is known to be censored (ξ_i = 1, δ_i = 0), then inverse probability weighting leads to synthetic data of the form Ŷ_i,W = {1 − 1/_n(Ỹ_i)}_n(Z_i)Ỹ_i/{1 − Ĝ_n(Ỹ_i)}, the sign of which is the opposite of Ỹ_i, and hence Y_i. Thus, this method creates synthetic data that are more variable than those based on either Y_R or Y_I, which can be seen by noting

Censoring rates of 20 and 40% were obtained by setting μ = 5.395 and 4.067, respectively. Values for θ₁ and θ₂ were selected by specifying the proportion of subjects with a missing censoring indicator at ỹ = 0 and by specifying the average proportion (over all ỹ) with a missing censoring indicator. We simulated data so that the average missingness rate was 20 or 40% and so that the missingness rate increased or decreased with ỹ. We chose θ₁ so that 1 − π(0) = {1 + exp(θ₁)}⁻¹ was 0.15 or 0.25 when the average missingness rate was 20%, and so that 1 − π(0) was 0.30 or 0.50 when the average missingness rate was 40%. We also examined situations in which none of the censoring indicators were missing, both for 20 and 40% censoring, as well as the case where there was no censoring at all.

We generated 10,000 Monte Carlo random samples of size n = 50, 100, 200 and 400 under each combination of censoring and missingness. Every data set was analyzed under the model: Y = α + βX + ε. For each configuration, we averaged over the 10,000 data sets to estimate the mean squared error (MSE), bias, and standard error (SE) associated with the slope estimators _R, _I, and _W. The results for the intercept estimators _R, _I, and _W were qualitatively the same, so they are not presented. When none of the censoring indicators were missing, we also calculated the unbounded Koul et al. (1981) estimator: where is defined in (1). We view _K as a reference for comparisons in our simulations.

We used the uniform kernel function if |u| ≤ 1 and W(u) = 0 otherwise, and the biweight kernel function if |u| ≤ 1 and K(u) = 0 otherwise. The bandwidths were h_n = b_n = n^−1/3max(ỹ). We estimated m(z) = pr(δ = 1|z) under the logistic model: log{m(z)/[1 − m(z)]} = γ₁ + γ₂ỹ + γ₃x. When the data on δ are completely (or quasi-completely) separated, the maximum likelihood estimate of γ = (γ₁, γ₂, γ₃) does not exist (Albert and Anderson 1984; Santner and Duffy 1986). We excluded such data sets and continued simulating until we obtained 10,000 samples. The proportion of samples excluded did not exceed 0.3%, except when the average missingness rate was 40%, the censoring rate was 20%, and n was 50, in which case less than 2.3% of the samples were excluded.

The MSE, bias, and SE for _K, _R, _I, and _W are presented in Tables 1, ​,2,2, and ​and3,3, respectively. The average jackknife estimate of SE is also shown in Table 3. Results are given only for situations where the average missingness rate increased with ỹ, as the results for the decreasing cases were virtually identical. The Koul et al. (1981) estimator (_k) is included as a benchmark when no censoring indicators are missing. The first row of each table corresponds to the special case of no censoring, where all of the estimators are identical because each synthetic response reduces to the observed response Y_i.

Table 1 shows that when no censoring indicators were missing, the MSEs for _R were less than or equal to those for _K, whereas the MSEs for _I and _W were slightly larger, at least for the 40% censoring case. For any given configuration, with or without missing censoring indicators, the MSE for _R never exceeded the MSE for _I or _W. As expected, the MSEs for _R, _I, and _W decreased (i.e., improved) as sample sizes increased, as censoring rates decreased, and as average missingness rates decreased (see Table 1).

The estimators also were evaluated with respect to bias and SE, the squares of which contribute equally to MSE. Table 2 gives the biases, all of which are negative when censoring is present; thus, each approach tended to underestimate the slope parameter. One interesting result is that applying our methods with missing censoring indicators usually produced less bias than applying the standard estimator of Koul et al. (1981) with complete censoring information. In fact, with a 40% censoring rate, the bias for each of our estimators (with missingness rates of either 20 or 40%) was always less than for _K with no missing censoring indicators. Among the proposed estimators, the bias was consistently smallest for the inverse probability weighted estimator, _W. Also, biases decreased as sample sizes increased and as censoring rates decreased, but were not affected as much by missingness rates.

Table 3 gives the SEs and their jackknife estimates. The SE patterns mimic those for the MSEs in Table 1 because bias and SE contribute equally to MSE and the SEs are larger (in absolute value) than the biases. As with the MSEs, _R had the smallest SEs and _W had the largest. The SEs decreased as n increased, as censoring decreased, and as missingness decreased. The jackknife estimates were always good for _W, but became too small for _I, and even smaller for _R, as n decreased, as censoring increased, and as missingness increased.

We have data on 276 glioblastoma patients who were over 50 years old when they entered the trial. By the end of the study, 59 of these patients progressed and were no longer ambulatory, 14 progressed and were still ambulatory, 166 did not progress, and 37 progressed but had an unknown ambulatory status. Therefore, with respect to the time to non-ambulatory progression, we have 59 uncensored observations (ξ = 1, δ = 1), 180 censored observations (ξ = 1, δ = 0), and 37 observations with a missing censoring indicator (ξ = 0). We have data on a covariate vector X = (X₁, X₂, X₃), where X₁ is fixed at 1, X₂ indicates the patient’s sex (0 for female, 1 for male), and X₃ is the patient’s age (in years) at study entry. Our data set includes information on 109 women and 167 men, with ages ranging from 51 to 74.

We applied our methods to the glioblastoma data to estimate the effects of sex and age on the time to non-ambulatory progression. As a first step, we fitted the linear regression model: E(Y|X = x) = β₁x₁ + β₂x₂ + β₃x₃, where Y is the log of the time to non-ambulatory progression. This initial model specifies parallel lines, with a separate intercept for each sex and a common slope in age. To help assess the appropriateness of the model, we added a sex-age interaction term (β₄x₂x₃), which allows each line to have its own sex-specific intercept and slope. Estimated regression coefficients and their standard errors, obtained under both models via all three proposed methods, are shown in Table 5. As a further means of assessing model suitability, we also fitted a separate quadratic model in age for each sex, as well as a model with quadratic-age curves that shared common age terms, but none of the quadratic terms were significant, so these results are not presented.

Based on 2-sided Wald tests, none of the analyses provided evidence of a sex effect, but to varying degrees they all suggested a detrimental age effect, with time to nonambulatory progression decreasing with age. The expanded model gave no hint of a sex-age interaction, regardless of which method was used, so we focus on the simpler model with only main effects for sex and age. Relative to the typical 0.05 level, the statistical significance of the age effect was marginal (p = 0.047) in the inverse probability weighted analysis, slightly stronger (p = 0.029) in the imputation analysis, and fairly high (p < 0.001) in the regression calibration analysis. Though excluding 13% of the data seems inefficient, if we ignore the 37 observations with a missing censoring indicator, the Koul et al. (1981) analysis gives similar results, with no evidence of a sex effect and a significant age effect (p = 0.022).

Thus, time to non-ambulatory progression does not differ significantly between men and women, but it appears to decrease (i.e., quality of life worsens) with age. Although evidence of this age effect is weak for the inverse probability weighted analysis, it is consistent across all three of the proposed analyses, as well as the Koul et al. (1981) analysis. The relative sizes of the estimated standard errors in Table 5 agree with our simulation results, where the SE for _R was smallest (and its jackknife estimates were too small), followed by the SE for _I, and finally by the SE for _W (with good jackknife estimates). This suggests that the test based on _W might be the least powerful, which further strengthens the conclusion that the observed age effect is real, as even the least powerful test was statistically significant.

Similar to Dikta (1998), our regression calibration method can be used even when all of the censoring indicators are observed. In that situation, the Koul et al. (1981) estimator (_K) transforms the observations (Ỹ_i, δ_i) into synthetic responses and then applies a modified least squares procedure to . Their approach is motivated by the fact that , but _K can perform poorly with small sample sizes. Under the MAR assumption, we have which implies that _K might be improved by replacing δ_i in with a regression estimator of m(Z_i) based on {(Ỹ_i, δ_i) : i = 1,…,n}. Our simulation results support this suggestion, but a detailed discussion is beyond the scope of this paper.

Along these same lines, we show in Sects. 3 and 4 that _I and _W have larger asymptotic variances than _R, and thus our imputation and inverse probability weighted estimators are asymptotically less efficient than our regression calibration estimator. The simulation studies confirm these results (see Tables 3 and ​and4).4). This suggests that, at least with respect to this measure of performance, the regression calibration method for handling missing censoring indicators is better than the other two methods, even though the other two methods can lead to asymptotically efficient estimators in some situations (see, e.g., Robins et al. 1994; Hahn 1998; Hirano et al. 2003; and Wang et al. 2004).

On the other hand, there are several good reasons to favor _W over _I or _R. The inverse probability weighted estimator is less biased than the other estimators (see Tables 2 and ​and4)4) and its jackknife estimates of SE were the most accurate (see Table 3). Also, the inverse probability weighted estimator enjoys the double-robustness property, which makes it less susceptible to problems caused by incorrectly specified models (see Table 4).

How should one balance the advantages and disadvantages of each approach? In general, we recommend the inverse probability weighted estimator. Not only is the inverse probability weighted estimator more robust to poor model choices for m(z), but relative to the other estimators, it is less biased and its jackknife variance estimates are more accurate. Although the regression calibration estimator has the smallest variance, its relatively large bias can produce misleading results in practice. In fact, our simulation showed that if a poor model was selected for m(z), the bias of _R could actually increase with the sample size. Currently, the variance estimator for _R does not adequately account for uncertainty in estimating m(z) and G(z). Also, the variance estimates are functions of the parameter estimates and _R has the largest bias; thus, even though the true variance of _R is smallest, the variance of _R also can be severely underestimated (see Table 3), making results appear too significant. Future efforts should focus on accounting for uncertainty in estimating m(z) and G(z), but until then, we cannot recommend using the regression calibration method in its present form.

Clearly, our estimators depend on the choice of bandwidths. However, each is a global functional and thus the rate n^1/2 asymptotic normality of the proposed estimators suggests that a proper choice of b_n in assumption (C.7) does not depend on the first order term of the mean squared error. Rather, it may depend only on second (or higher) order terms. This implies that the selection of b_n may be fairly flexible as long as the bandwidths satisfy (C.7).