Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Lifetime Data Anal. Author manuscript; available in PMC 2012 April 1.
Published in final edited form as:
PMCID: PMC3020262

Linear regression analysis of survival data with missing censoring indicators


Linear regression analysis has been studied extensively in a random censorship setting, but typically all of the censoring indicators are assumed to be observed. In this paper, we develop synthetic data methods for estimating regression parameters in a linear model when some censoring indicators are missing. We define estimators based on regression calibration, imputation, and inverse probability weighting techniques, and we prove all three estimators are asymptotically normal. The finite-sample performance of each estimator is evaluated via simulation. We illustrate our methods by assessing the effects of sex and age on the time to non-ambulatory progression for patients in a brain cancer clinical trial.

Keywords: Asymptotic normality, Censoring indicator, Imputation, Inverse probability weighting, Least squares, Missing at random, Regression calibration

1 Introduction

Linear regression is a popular statistical tool that has been used successfully in many areas. In survival analysis, a log-transformation of the response variable converts a conventional linear model to an accelerated failure time model, which is an appealing alternative to the Cox (1972) proportional hazards model because of its direct interpretation (cf. Reid 1994). Linear regression has been studied extensively for the analysis of randomly right censored data, and various procedures have been developed; see, for example, Miller (1976), Buckley and James (1979), Koul et al. (1981), Leurgans (1987), Ritov (1990), Tsiatis (1990), Lai and Ying (1991); Ying (1993), Jin et al. (2003), and Li and Wang (2003). All of these methods assume that every censoring indicator is observed.

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

Yi=Xi[top top]β+εi,

where Yi is a response variable, Xi 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 {X1,…,Xn}, the errors {ε1,…,εn} are independent and identically distributed (i.i.d.) with mean zero. Let {C1,…,Cn} be i.i.d. censoring variables, where Ci can randomly censor Yi on the right. The observed response is i = min(Yi, Ci), the censoring indicator is δi = I (YiCi), 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 Yi*, where Gn(c) is a product-limit type estimator of the censoring distribution G(c) = pr(Cc). Koul et al. estimated β by applying a modified least squares procedure to (Yi*,Xi), but their estimator does not handle missing censoring indicators. Also, they deflate the synthetic response Yi* to 0 when the true response Yi exceeds Ci, and they inflate Yi* beyond Yi when Yi 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 mi, say piδi + (1 − pi)mi, which allows statistical inference about β when some censoring indicators are missing.

The article is organized as follows. Section 2 develops a synthetic data approach based on regression calibration, which sets pi [equivalent] 0. Sections 3 and 4 define imputation and inverse probability weighted estimators, which set pi = ξi and pi = ξi/Êi|i), respectively. As δi is unknown if ξi = 0, each of the three formulas for pi is zero in that situation. All three least squares estimators of β are shown to be asymptotically normal. Section 5 reports the results of a simulation study conducted to evaluate the finite sample behavior of the three proposed estimators, and Sect. 6 illustrates their use by applying them to the glioblastoma data. Finally, Sect. 7 discusses our estimators and the underlying assumptions, and proofs of the main results are given in the Appendix.

2 Least squares regression calibration estimator

If we define Zi[top top]=(Y˜i,Xi[top top]) and m(zi) = Ei|Zi = zi), then one natural generalization of Yi* replaces δi by its conditional expectation to obtain Yi,R = m(Zi)i/{1 − G (i)}. Averaging over δi and invoking the usual conditional expectation arguments yield:

E(Yi,R|Xi)=(...)=E[E{Y˜iδi1G(Y˜i)|Y˜i,Xi}|Xi]=(...)=E(Yi|Xi)=Xi[top top]β.

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

β˜R=(i=1nXiXi[top top])1i=1nXiYi,R.

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

For example, consider a parametric model for m(z) of the form m0(z, θ), where m0(·, ·) is a known continuous function and θ is a vector of unknown parameters. The parameter vector θ can be estimated by [theta w/ hat]n, the maximizer of the likelihood


and thus m(z) can be estimated by mn(z) = m0(z, [theta w/ hat]n). As for estimating G(·), we begin by defining π() = E (ξ| = ), which can be estimated nonparametrically by


where W(·) is a kernel function and bn 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 hn is a bandwidth sequence. Finally, similar to Dikta (1998) and Wang and Ng (2008), we define the following estimator of G():


where Ri denotes the rank of i (i = 1,…,n). Now, analogous to betaR, we can define the following least squares regression calibration estimator of β:

β^R=(i=1nXiXi[top top])1i=1nXiY^i,R,

where Ŷi,R is Yi,R with m(·) and G(·) replaced by mn(·) and Ĝn(·), respectively.

Define σr,s=E[π(Y˜)m0(Z,θ){1m0(Z,θ)}[partial differential]m0(Z,θ)[partial differential]θr[partial differential]m0(Z,θ)[partial differential]θs] for 1 ≤ r, sd + 1 and assume:

  • (C.1) m0(z, θ) has bounded first-order derivatives on θ.
  • (C.2) σr,s < ∞ for 1 ≤ r, sd + 1 and matrix A(θ) = (σr,s) is positive definite.
  • (C.3) E {‖X22/[H with macron]2()} < ∞, where H(t) = pr(t) and [H with macron](t) = 1 − H (t).
  • (C.4) E {‖X2} < ∞.
  • (C.5) π(·) has a bounded first-order derivative.
  • (C.6) W(·) is of order 1 with bounded support.
  • (C.7) nbn → ∞ and nbn20.

The following theorem addresses the asymptotic normality of [beta]R.

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


where VR = Σ−1RΣ−1, Σ = E(XX[top top]), andR is defined in (A.15) of the Appendix.

The structure ofR is complex; thus, estimating VR by replacing the unknowns inR with appropriate component estimates is difficult to implement in practice. Alternatively, we estimate VR by nVJ, where VJ is a jackknife estimator of the asymptotic variance of [beta]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:

VJ^=(i=1nXiXi[top top])1[i=1n(Y^i,RXi[top top]β^R)2XiXi[top top]](i=1nXiXi[top top])1.

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


Proofs of Theorem 2.1 and Theorem 2.2 are given in the Appendix.

3 Least squares imputation estimator

Imputation is a popular method for dealing with missing data, and this section describes such an approach. When some censoring indicators are missing, some synthetic data in (1) are not defined. However, in the expression for Yi*, we can replace Gn(·) by Ĝn(·) and impute δi with mn(Zi) if δi is missing (ξi = 0), which yields the imputed synthetic data

Y^i,I={ξiδi+(1ξi)m^n(Zi)}Y˜i1G^n(Y˜i) i=1,(...),n.

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

β^I=(i=1nXiXi[top top])1i=1nXiY^i,I.

Under the assumed missingness mechanism, [beta]I also can be motivated by the relationship

E(YI|X)=E([E(ξ|Y˜)δ+{1E(ξ|Y˜)}m(Z)]Y˜21G(Y˜)|X)=(...)=E(YR|X)=X[top top]β,

where YI = {ξδ + (1 − ξ)m(Z)}/{1 − G()} and YR = m(Z)/{1 − G()}.

The following theorem addresses the asymptotic normality of [beta]I.

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


where VI = Σ−1IΣ−1, ΩI=ΩR+ΩI*, and ΩI* is defined in (A.20) of the Appendix.

Similar to Sect. 2, we estimate VI by the jackknife methods described earlier. Clearly [beta]I has larger asymptotic variance than [beta]R, and hence [beta]I is asymptotically less efficient than [beta]R. Intuitively, one reason for this may be that the imputed synthetic observation Ŷi,I is zero if Yi is known to be censored (ξi = 1, δi = 0). We see that the imputation approach creates synthetic data with greater variability than the regression calibration method by noting


4 Least squares inverse probability weighted estimator

Inverse probability weighting is another popular approach for dealing with missing data. If we define YW = [{ξ/π()}δ + {1 − ξ/π()}m(Z)]/{1 − G()}, then under the assumed missingness mechanism, conditional expectations similar to those in Sects. 2 and 3 give

E(YW|X)=E([{E(ξ)/π(Y˜)}δ+{1E(ξ)/π(Y˜)}m(Z)]Y˜1G(Y˜)|X)=(...)=E(YR|X)=X[top top]β.

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

Y^i,W=[{ξiπ^n(Y˜i)}δi+{1ξiπ^n(Y˜i)}m^n(Zi)]Y˜i1G^n(Y˜i) i=1,(...),n.

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

β^W=(i=1nXiXi[top top])1i=1nXiY^i,W.

The following theorem addresses the asymptotic normality of [beta]W.

Theorem 4.1 Under assumptions (C.1) through (C.7), we have


where VW = Σ−1WΣ−1, ΩW=ΩR+ΩW*, and ΩW* is defined in (A.26) of the Appendix.

As described earlier, VW can be estimated by jackknife methods. We see [beta]W has a larger asymptotic variance than [beta]I, and hence [beta]R, by noting α[top top](ΩW*ΩI*)α>0 for any d-vector α. One explanation is that if Yi is known to be censored (ξi = 1, δi = 0), then inverse probability weighting leads to synthetic data of the form Ŷi,W = {1 − 1/[pi]n(i)}mn(Zi)i/{1 − Ĝn(i)}, the sign of which is the opposite of i, and hence Yi. Thus, this method creates synthetic data that are more variable than those based on either YR or YI, which can be seen by noting


One advantage of this approach is that if either m(·) or π(·) is modeled correctly, the estimator [beta]W will be consistent, which is the so-called “double-robustness” property (see, e.g., Scharfstein et al. 1999; Lunceford and Davidian 2004; and Wang et al. 2004). Note that we estimate π(·) nonparametrically, and thus [beta]W enjoys this robustness even if the model for m(·) is not specified correctly. This is an appealing property of [beta]W, but must be weighed against its larger asymptotic variance.

5 Simulation study

The finite-sample properties of the proposed estimators were evaluated via Monte Carlo simulation, with emphasis on the effects of sample size, percentage of observations censored, and percentage of censoring indicators missing. Response data were generated from the linear model: Y = α + βX + ε, with α = 3 and β = 0.5. The random variables ε, X, and C were independently generated from N(0,1), U(0,2), and N(μ, 4) distributions, respectively. For each subject, the observed response was = min(Y, C) and the censoring indicator was δ = I (YC). Conditional on = , the probability that the censoring indicator was missing, 1 − π(), was determined by the logistic model: log{π()/|1 − π()]} = θ1 + θ2.

Censoring rates of 20 and 40% were obtained by setting μ = 5.395 and 4.067, respectively. Values for θ1 and θ2 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 θ1 so that 1 − π(0) = {1 + exp1)}−1 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 [beta]R, [beta]I, and [beta]W. The results for the intercept estimators [alpha]R, [alpha]I, and [alpha]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:

β^K=(i=1nXiXi[top top])1i=1nXiYi*,

where Yi* is defined in (1). We view [beta]K as a reference for comparisons in our simulations.

We used the uniform kernel function W(u)=12 if |u| ≤ 1 and W(u) = 0 otherwise, and the biweight kernel function K(u)=1516(12u2+u4) if |u| ≤ 1 and K(u) = 0 otherwise. The bandwidths were hn = bn = n−1/3max(). We estimated m(z) = pr(δ = 1|z) under the logistic model: log{m(z)/[1 − m(z)]} = γ1 + γ2 + γ3x. When the data on δ are completely (or quasi-completely) separated, the maximum likelihood estimate of γ = (γ1, γ2, γ3) 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 [beta]K, [beta]R, [beta]I, and [beta]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 ([beta]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 (Yi*,Y^i,R,Y^i,I,Y^i,W) reduces to the observed response Yi.

Table 1
Mean squared errors for the regression calibration ([beta]R), imputation ([beta]I), and inverse probability weighted ([beta]W) estimators of β by censoring rate (CR), missingness rate (MR), and sample size (n) ...
Table 2
Biases for the regression calibration ([beta]R) imputation ([beta]I), and inverse probability weighted ([beta]W) estimators of β by censoring rate (CR), missingness rate (MR), and sample size (n) based on 10,000 ...
Table 3
Standard errors (and mean jackknife estimates) for the regression calibration ([beta]R), imputation ([beta]I), and inverse probability weighted ([beta]W) estimators of β by censoring rate (CR), missingness rate ...

Table 1 shows that when no censoring indicators were missing, the MSEs for [beta]R were less than or equal to those for [beta]K, whereas the MSEs for [beta]I and [beta]W were slightly larger, at least for the 40% censoring case. For any given configuration, with or without missing censoring indicators, the MSE for [beta]R never exceeded the MSE for [beta]I or [beta]W. As expected, the MSEs for [beta]R, [beta]I, and [beta]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 [beta]K with no missing censoring indicators. Among the proposed estimators, the bias was consistently smallest for the inverse probability weighted estimator, [beta]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, [beta]R had the smallest SEs and [beta]W had the largest. The SEs decreased as n increased, as censoring decreased, and as missingness decreased. The jackknife estimates were always good for [beta]W, but became too small for [beta]I, and even smaller for [beta]R, as n decreased, as censoring increased, and as missingness increased.

In order to assess robustness, we analyzed the same simulated data with a “poor” model choice for m(z). Rather than using log{m(z)/[1 − m(z)]} = γ1 + γ2 + γ3x to model m(z) as a logistic function of and x, we set γ2 = γ3 = 0 and modeled m(z) as a constant, which gives m^(z)[equivalent]m^=i=1nξiδi/i=1nξi. Table 4 displays the biases and SEs obtained for this poor model choice. The results for no missing censoring indicators (M R = 0%) are the same as in Tables 2 and and3,3, except for [beta]R, as only [beta]R depends on m(z) in this situation. The SEs always decrease as sample sizes increase, as does the bias of [beta]W, but the biases of [beta]R and [beta]I actually increase with the sample size in many cases. In fact, the bias of [beta]R is now positive and increases with n in every situation, whereas it was always negative and approached 0 as n increased when using the better model for m(z). The bias of [beta]I increases with the sample size, which makes the negative biases smaller in absolute value, but once the biases become positive, they grow in absolute value as n increases. In contrast, the bias of [beta]W for the poor model choice is nearly identical to the bias of [beta]W when using the better model. As predicted theoretically, the double-robustness property of [beta]W apparently keeps its performance from degrading if a poor choice is made when modeling m(z).

Table 4
Biases (and standard errors) for the regression calibration ([beta]R), imputation ([beta]I), and inverse probability weighted ([beta]W) estimators of β for a poor choice of m(z) [i.e., m(z) [equivalent] m] by censoring ...

6 Example: analysis of glioblastoma data

We illustrate our methods with the glioblastoma data introduced earlier. These data are from a brain cancer clinical trial. As a measure of post-treatment quality of life, we focus on the number of days until a patient declines to a state in which he has lost his mobility and his cancer has returned. Thus, the endpoint of interest is the time (in days) to non-ambulatory progression, and we want to relate that outcome to the patient’s sex and age at study entry.

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 = (X1, X2, X3)[top top], where X1 is fixed at 1, X2 indicates the patient’s sex (0 for female, 1 for male), and X3 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) = β1x1 + β2x2 + β3x3, 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 (β4x2x3), 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.

Table 5
Estimated regression coefficients (and jackknife estimates of their standard errors) under two models for the glioblastoma data from the brain cancer clinical trial (n = 276)

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 [beta]R was smallest (and its jackknife estimates were too small), followed by the SE for [beta]I, and finally by the SE for [beta]W (with good jackknife estimates). This suggests that the test based on [beta]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.

In this example, expressing the failure time of interest as a function of two component times reveals a potential source of bias. That is, the time to non-ambulatory progression can be viewed as the maximum of the time to loss of mobility and the time to cancer progression. We treat patients who have progressed but are still ambulatory at the end of the study as contributing censored failure times, but this censoring is likely informative and may introduce a bias, as ambulatory patients who have already progressed should be at greater risk of failing than ambulatory patients who have not yet progressed. Similarly, patients in remission who have already lost their mobility should be at greater risk of failing than patients in remission who are still ambulatory. This problem is not unique to our approach, however, and affects a wide class of existing survival techniques, including something as simple as using a Kaplan–Meier curve to estimate the distribution of time to non-ambulatory progression. Future efforts should be directed at developing methods that properly account for this type of informative censoring, such as a multivariate procedure, although such developments are beyond the scope of this article.

7 Discussion

We propose three estimators for the regression coefficients in a linear least squares analysis of censored survival data when some censoring indicators are missing. Our methods generalize the synthetic data methods of Koul et al. (1981), which do not allow for missing censoring indicators. The proposed estimators modify the synthetic responses of Koul et al. so that each missing censoring indicator (δ) is replaced by its estimated conditional expectation. Our estimators differ from each other only with respect to how they handle the known censoring indicators. The regression calibration estimator replaces each known δ by its estimated mean, the imputation estimator uses the known δ, and the inverse probability weighted estimator substitutes a weighted average of the known δ and its estimated mean.

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 ([beta]K) transforms the observations (i, δi) into synthetic responses (Yi*) and then applies a modified least squares procedure to (Yi*,Xi). Their approach is motivated by the fact that E[δiY˜i{1G(Y˜i)}1|Xi]=Xi[top top]β, but [beta]K can perform poorly with small sample sizes. Under the MAR assumption, we have


which implies that [beta]K might be improved by replacing δi in Yi* with a regression estimator of m(Zi) 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 [beta]I and [beta]W have larger asymptotic variances than [beta]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 [beta]W over [beta]I or [beta]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 [beta]R could actually increase with the sample size. Currently, the variance estimator for [beta]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 [beta]R has the largest bias; thus, even though the true variance of [beta]R is smallest, the variance of [beta]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 [beta] is a global functional and thus the rate n1/2 asymptotic normality of the proposed estimators suggests that a proper choice of bn 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 bn may be fairly flexible as long as the bandwidths satisfy (C.7).

We assume the censoring variables (C1,…,Cn) are i.i.d., which is the same assumption made by Koul et al. (1981). This form of censoring is more restrictive than the standard noninformative censoring assumed in most survival analyses, as the latter does not require the censoring distribution G(c) to be estimated. Our methods, however, as well as those developed by Koul et al., rely on synthetic data that are functions of an estimate of G(c). One avenue of future research might involve extending our methods by expressing G(c) as a function of the covariates. Generally, this requires either a parametric model, the choice of which must be justified, or else a nonparametric approach, which could suffer in high-dimensional situations if the data were too sparse. Alternatively, we might view G(c) as the marginal distribution of C, which “averages” over the covariates.


Qihua Wang’s research was supported by the National Science Fund for Distinguished Young Scholars in China (10725106), the National Natural Science Foundation of China (10671198), the National Science Fund for Creative Research Groups in China, the research Grants Council of Hong Kong (HKU 7050/06P), and the grant from Key Lab of Random Complex Structures and Data Science, CAS. Gregg Dinse’s research was supported by the Intramural Research Program of the NIH, National Institute of Environmental Health Sciences (Z01-ES-045007-13). The authors thank Dr. Shyamal Peddada and Dr. David Dunson, as well as the editors and referees, for their many helpful comments.

Appendix: Proofs of theorems

Define Grad[m0(z,θ)]=([partial differential]m0(z,θ)[partial differential]θ1,[partial differential]m0(z,θ)[partial differential]θ2,,[partial differential]m0(z,θ)[partial differential]θd,[partial differential]m0(z,θ)[partial differential]θd+1)[top top] and set A(θ)=(σr,s) for 1≤r, sd+1, where σr,s=E[π(Y˜)m0(Z,θ){1m0(Z,θ)}[partial differential]m0(Z,θ)[partial differential]θr[partial differential]m0(Z,θ)[partial differential]θs]. Furthermore, let α(z1, z2) = Grad[m0(z1, θ)][top top] A−1(θ) Grad[m0(z2, θ)] and assume the following conditions:

  • (C.1) m0(z, θ) has bounded first-order derivatives on θ.
  • (C.2) σr,s < ∞ for 1 ≤ r, sd + 1 and the matrix A(θ) = (σr,s) is positive definite.
  • (C.3) E {‖X22/[H with macron]2()} < ∞, where H(t) = pr(t) and [H with macron](t) = 1 − H(t).
  • (C.4) E {‖X2} < ∞.
  • (C.5) π(·) has a bounded first-order derivative.
  • (C.6) W(·) is of order 1 with bounded support.
  • (C.7) nbn → ∞ and nbn20.

Proof of Theorem 2.1 For Σn   =   1ni=1nXiXi[top top] and An   =   n1/2i=1nXi(Y^i,RXi[top top]β), we can write


By adding and subtracting identical terms, we can rewrite An as follows:

An=n1/2i=1nXi(Yi,RXi[top top]β)+n1/2i=1nXiY˜i{m0(Zi,θ^n)m0(Zi,θ)}1G(Y˜i)+n1/2i=1nXi{Y˜im0(Zi,θ)1G(Y˜i)}G^n(Y˜i)G(Y˜i)1G(Y˜i)+op(1):=Tn1+Tn2+Tn3+op(1).

Define μ(z)=E{XY˜1G(Y˜)α(Z,Zj)|Zj=z}. Under conditions (C.1) and (C.2), we have


Set H0(t) = pr(i > t, δi = 0) and define

ψ(Y˜i,δi,ξi;t)=[{ξiπ(Y˜i)}{δiu(Y˜i)}π(Y˜i){1H(Y˜i)}I(Y˜it)+0Y˜i[logical or operator]tdH˜0(S){1H(s)}2+I(Y˜it,δi=0)1H(Y˜i)],

where u(·) is given in Sect. 2. Similar to Wang and Ng (2008), we have


The third term in (A.2) can be rewritten as




and set Un = n−3/2i<j h(Xi, i, δi, ξi; Xj, j, Sj, ξj), which leads to


Combining (A.2), (A.3), and (A.6) gives


where Tn1 is defined in (A.2) and Mn is the main component of Tn2 in (A.3).

By the central limit theorem, we have


and under the assumed missingness mechanism, we have


where Ω1 = E{XX[top top](YR − X[top top]β)2} and Ω2=E[μ(Y˜)μ[top top](Y˜)π(Y˜)m0(Z,θ){1m0(Z,θ)}]. For any d-vector of constants, say α, let hα(·) = α[top top]h(·), Unα = α[top top]Un, and gα(·) = α[top top]g(·), where


Under the assumed missingness mechanism, we have

E{ψ(Y˜i,δi,ξi;Y˜j)|Y˜j}=0  forij,

from which it is straightforward to verify that E{hα(X1, 1, δ1, ξ1; X2, 2, δ2, ξ2)} = 0. As a result of the relationship

E{ψ2(Y˜i,δi,ξi;Y˜j)|Y˜j}={1G(Y˜j)}20Y˜j{H¯(s)}2dH˜0(s){1G(Y˜j)}2H¯2(Y˜j)   for  ij

and condition (C.3), it follows that E(hα2)<. By the central limit theorem for U-statistics, we have

UnαLN(0,α[top top]Ω3α),


Ω3=E{g(X1,Y˜1,δ1,ξ1)g[top top](X1,Y˜1,δ1,ξ1)}

and α[top top]Ω3α=E(gα2). This result implies


As for covariances, under the assumed missingness mechanism, we have

cov(Tn1,Mn)=E[Xi{Y˜im0(Zi)1G(Zi)Xi[top top]β}μ[top top](Y˜i)π(Y˜i){δim0(Y˜i)}m0(Zi,θ){1m0(Zi,θ)}]=0

Note that E{Xk(Yk,RXk[top top]β)hij[top top]}=0 for ki, kj, i ≠ j and E{Xi(Yi,RXi[top top]β)hij}=E{Xj(Yj,RXj[top top]β)hij}, where hij = h(Xi, i, δi, ξi; Xj, j, δj, ξj).

Similarly, we have

cov(Tn1,Un)=(11n)E{E[X1(Y1,RX1[top top]β)ψ(Y˜1,δ1,ξ1;Y˜2)|Y˜2]X2[top top]Y˜2m0(Z2,θ)1G(Y˜2)}E{E[X1(Y1,RX1[top top]β)ψ(Y˜1,δ1,ξ1;Y˜2)|Y˜2]X2[top top]Y˜2m0(Z2,θ)1G(Y˜2)}=E{X1(Y1,RX1[top top]β)X2[top top]Y˜2m0(Z2,θ)1G(Y˜2)[0Y˜2[logical or operator]Y˜1dH˜0(s){1H(s)}2+I(Y˜1Y˜2,δ1=0)1H(Y˜1)]}:=Ω1,3.


cov(Mn,Un)=(11n)E[m0(Z2,θ)Y˜2X2μ[top top](Y1){1G(Y˜2)}{1H(Y˜1)}I(Y˜1Y˜2)]E[m0(Z2,θ)Y˜2X2μ[top top](Y˜1){1G(Y˜2)}{1H(Y˜1)}I(Y˜1Y˜2)]:=Ω2,3.

Together Eqs. (A.7)(A.9) and (A.10)(A.13) prove




By the law of large numbers and (C.4), it follows that n1Σi=1nXiXi[top top]pΣ. This result, together with (A.1) and (A.14), prove Theorem 2.1.

Proof of Theorem 2.2 Asymptotic representations similar to (A.7) can be obtained for Σi=1n(Y^i,RXi[top top]β^R)2XiXi[top top] for i = 1,…,n. Then, by applying these representations to (2) and using some algebra, Theorem 2.2 can be proved.

Proof of Theorem 3.1 By standard arguments, we have

 n1/2i=1nXi(Y^i,IXi[top top]β)=n1/2i=1nXi[{δiξi+(1ξi)m0(Zi,θ)}Y˜i1G(Y˜i)Xi[top top]β]+n1/2i=1nXi(1ξi){m0(Zi,θ^n)m0(Zi,θ)}Y˜i1G(Y˜i)+n1/2i=1nXiY˜{δiξi+(1ξi)m0(Zi,θ)}1G(Y˜i){G^n(Y˜i)G(Y˜i)1G(Y˜i)}+op(1)=Tn1+Tn2+Tn3+Tn4+Tn5+Tn6+op(1),

where Tn1, Tn2 and Tn3 are defined in (A.2), and Tn4, Tn5 and Tn6 are defined as follows


By (A.2) and (A.7), we have


Let μ˜(Zj)=E{XY˜π(Y˜)1G(Y˜)α(Z,Zj)|Zj}. Under conditions (C.1) and (C.2), using arguments similar to those leading to (A.3)(A.5), we get


Recalling the definition of Tn4, it follows from (A.18) that


By the central limit theorem, we have



ΩI*=E[L(X,Y˜)L[top top](X,Y˜)π(Y˜)m0(Z,θ){1m0(Z,θ)}]

and L(X,Y˜)=XY˜1G(Y˜)μ˜(Y˜)m0(Z,θ){1m0(Z,θ)}. Similar to the way (A.5) was derived, we obtain


Note that


From the fact that E{ψ(j, δj, ξj; i)|i} = 0, it follows that


Equations (A.21) and (A.22) together prove that


Under the assumed missingness mechanism, a little algebra shows cov(Tn1, Tn45) = 0, cov(Un, Tn45) = 0 and cov(Mn, Tn45) = 0. This result, Eqs. (A.14), (A.16), (A.17), (A.19), and (A.23), and the fact that 1nΣi=1nXiXi[top top]pΣ together prove Theorem 3.1.

Proof of Theorem 4.1 Similar to the derivation of (A.16), by (C.5), (C.6) and (C.7), we have

n1/2i=1n(Y^i,WXi[top top]β)=n1/2Xi(Y^i,RXi[top top]β)+n1/2i=1nXiξiπ^n(Y˜i){δim0(Y˜i,θ)}Y˜i1G^n(Y˜i)+n1/2i=1nXiξiπ^n(Y˜i){m0(Zi,θ)m0(Zi,θn^)}Y˜i1G^n(Y˜i)=Tn1+Mn+Un+Sn+Rn+op(1),

where Tn1, Mn and Un are defined in (A.7), and Sn and Rn are given by


Similar to the way (A3) was derived, we can show that


Hence, it follows from the central limit theorem that



ΩW*=E[L˜(X,Y˜)L˜[top top](X,Y˜)m0(Z,θ){1m0(Z,θ)}π(Y˜)]

and L˜(X,Y˜)=XY˜1G(Y˜)π(Y˜)μ(Y˜)m0(Z,θ){1m0(Z,θ)}.

It can be shown that cov(Tn1 + Mn + Un, Sn + Rn) = 0. This result, together with Eqs. (A.7), (A.14), (A.24), and (A.25), and the fact that 1nΣi=1nXiXi[top top]pΣ, prove Theorem 4.1.

Contributor Information

Qihua Wang, Department of Mathematics and Statistics, Yunnan University, Kunming 650091, China. Academy of Mathematics and Systems Science, Chinese Academy of Science, Beijing 100190, China.

Gregg E. Dinse, Biostatistics Branch, National Institute of Environmental Health Sciences, Research Triangle Park, North Carolina 27709, USA.


  • Albert A, Anderson JA. On the existence of maximum likelihood estimates in logistic regression models. Biometrika. 1984;71:1–10.
  • Buckley J, James I. Linear regression with censored data. Biometrika. 1979;66:429–436.
  • Cox DR. Regression models and life tables (with discussion) J R Stat Soc B. 1972;34:187–220.
  • Dewanji A. A note on a test for competing risks with missing failure type. Biometrika. 1992;79:855–857.
  • Dikta G. On semiparametric random censorship models. J Stat Plan Inference. 1998;66:253–279.
  • Dinse GE. Nonparametric estimation for partially-complete time and type of failure data. Biometrics. 1982;38:417–431. [PubMed]
  • Gao G, Tsiatis AA. Semiparametric estimators for the regression coefficients in the linear transformation competing risks model with missing cause of failure. Biometrika. 2005;92:875–891.
  • Goetghebeur EJ, Ryan L. A modified logrank test for competing risks with missing failure type. Biometrika. 1990;77:207–211.
  • Goetghebeur EJ, Ryan L. Analysis of competing risks survival data when some failure types are missing. Biometrika. 1995;82:821–833.
  • Hahn J. On the role of the propensity score in efficient semiparametric estimation of average treatment effects. Econometrica. 1998;66:315–331.
  • Hirano K, Imbens GW, Ridder G. Efficient estimation of average treatment effects using the estimated propensity score. Econometrica. 2003;71:1161–1189.
  • Horvitz DG, Thompson DJ. A generalization of sampling without replacement from a finite universe. J Am Stat Assoc. 1952;47:663–685.
  • Jin Z, Lin DY, Wei LJ, Ying Z. Rank-based inference for the accelerated failure time model. Bio-metrika. 2003;90:341–353.
  • Koul H, Susarla V, van Ryzin J. Regression analysis with randomly right-censored data. Ann Stat. 1981;9:1276–1288.
  • Lai TL, Ying Z. Rank regression methods for left-truncated and right-censored data. Ann Stat. 1991;19:531–556.
  • Leurgans S. Linear models, random censoring and synthetic data. Biometrika. 1987;74:301–309.
  • Li G, Wang QH. Empirical likelihood regression analysis for right censored data. Stat Sinica. 2003;13:51–68.
  • Lo S-H. Estimating a survival function with incomplete cause-of-death data. J Multivar Anal. 1991;39:217–235.
  • Lu W, Liang Y. Analysis of competing risks data with missing cause of failure under additive hazards model. Stat Sinica. 2008;18:219–234.
  • Lu K, Tsiatis AA. Multiple imputation methods for estimating regression coefficients in the competing risks model with missing cause of failure. Biometrics. 2001;57:1191–1197. [PubMed]
  • Lunceford JK, Davidian M. Stratification and weighting via the propensity score in estimation of causal treatment effects: a comparative study. Stat Med. 2004;23:2937–2960. [PubMed]
  • McKeague IW, Subramanian S. Product-limit estimators and Cox regression with missing censoring information. Scand J Stat. 1998;25:589–601.
  • Miller RG. Least squares regression with censored data. Biometrika. 1976;63:449–464.
  • Peddada SD, Patwardhan G. Jackknife variance estimators in linear models. Biometrika. 1992;79:654–657.
  • Reid N. A conversation with Sir David Cox. Stat Sci. 1994;9:439–455.
  • Ritov Y. Estimation in a linear regression model with censored data. Ann Stat. 1990;18:303–328.
  • Robins JM, Rotnitzky A, Zhao LP. Estimation of regression coefficients when some regressors are not always observed. J Am Stat Assoc. 1994;89:846–866.
  • Santner TJ, Duffy DE. A note on A. Albert and J. A. Anderson’s conditions for the existence of maximum likelihood estimates in logistic regression models. Biometrika. 1986;73:755–758.
  • Scharfstein DO, Rotnitzky A, Robins JM. Adjusting for nonignorable drop-out using semiparametric nonresponse models (with discussion) J Am Stat Assoc. 1999;94:1096–1146.
  • Subramanian S. Asymptotically efficient estimation of a survival function in the missing censoring indicator model. Nonparametr Stat. 2004;16:797–817.
  • Subramanian S. Survival analysis for the missing censoring indicator model using kernel density estimation techniques. Stat Methodol. 2006;3:125–136. [PMC free article] [PubMed]
  • Tsiatis AA. Estimating regression parameters using linear rank tests for censored data. Ann Stat. 1990;18:354–372.
  • Tsiatis AA, Davidian M, McNeney B. Multiple imputation methods for testing treatment differences in survival distributions with missing cause of failure. Biometrika. 2002;89:238–244.
  • Wang QH, Ng K. Asymptotically efficient product-limit estimators with censoring indicators missing at random. Stat Sinica. 2008;18:749–768.
  • Wang QH, Linton O, Härdle W. Semiparametric regression analysis with missing response at random. J Am Stat Assoc. 2004;99:334–345.
  • Ying Z. A large sample study of rank estimation for censored regression data. Ann Stat. 1993;21:76–99.
  • Zhou X, Sun L. Additive hazards regression with missing censoring information. Stat Sinica. 2003;13:1237–1257.