An empirical likelihood method for semiparametric linear regression with right censored data.

This paper develops a new empirical likelihood method for semiparametric linear regression with a completely unknown error distribution and right censored survival data. The method is based on the Buckley-James (1979) estimating equation. It inherits some appealing properties of the complete data empirical likelihood method. For example, it does not require variance estimation which is problematic for the Buckley-James estimator. We also extend our method to incorporate auxiliary information. We compare our method with the synthetic data empirical likelihood of Li and Wang (2003) using simulations. We also illustrate our method using Stanford heart transplantation data.

1. Introduction

Suppose that one observes right censored regression data consisting of n i.i.d. triples (X , Z , δ ) = (X , Y ∧C , I[Y ≤ C ]), i = 1,…, n, where for subject i, Y is a known monotone transformation of the survival time of interest, C is the corresponding censoring time, and X = (X ,…, X ) is a vector of p covariates. We consider the problem of making inferences for the slope parameter of the semiparametric linear regression model as follows: where β = (β 1,…, β ) is a p × 1 vector of unknown regression coefficients and ϵ 's are independent and identically distributed random errors with an unknown distribution F. Assume that ϵ is independent of (X , C ) and that conditional on X , Y and C are independent, i = 1,…, n. Because the error distribution F is completely unknown, we assume no intercept term in model (1) in order for F to be identifiable. Assume further that the covariance matrix of X is positive definite. Model (1) provides a useful alternative to the popular Cox [1] model when the proportional hazards assumption does not hold. Furthermore, it makes no parametric assumption on the error distribution and is, thus, more flexible than parametric accelerated failure time models (cf. Andersen et al. [2]).

Buckley and James [3] extended the least squares method to estimate the regression coefficient β in model (1) with right censored data. The Buckley-James estimator can be calculated using an iterative algorithm. Its consistency and asymptotic normality have been established by Lai and Ying [4], Ritov [5], Tsiatis [6], Ying [7], and others. However, its asymptotic variance involves the unknown hazard function of ϵ, and its derivatives whose nonparametric estimations are problematic (cf. Ritov [5]). To overcome this problem, several approaches have been studied in the literature. Jin et al. [3, 8], Wei et al. [9], and Lin and Geyer [10] considered rank-based inferences. Koul et al. [11] introduced a synthetic data approach which was further investigated by Zhou [12], Lai et al. [13], and others. Let Y = δ Y /G(Y −), where G is the survival function of the censoring time. It can be shown that E(Y | X ) = E(Y | X ) if C is independent of X and Y. Koul et al. [11] proposed to estimate the regression parameters by regressing on X , where is the Kaplan-Meier [14] estimate of G and is referred to as a synthetic variable. Koul et al. [11] showed that the synthetic data estimator is asymptotically normal and its variance is simple to estimate. However, it requires a strong assumption that the censoring time is independent of both the survival time and the covariates. It can also have poor small sample performance in the presence of heavy censoring (cf. Li and Wang [15]). Jing and Qin [16] and Li and Wang [15] independently developed empirical likelihood-based inference for β using synthetic data. Their methods provide substantial improvement over the normal theory method of Koul et al. [11], but are still not very satisfactory for small samples (cf. Li and Wang [15]). The synthetic data approach makes an independence assumption between C and (Y, X). Li and Lu [17] showed that it can yield seriously biased parameter estimate if C depends on X. Zhou and Li [18] developed a censored data EL method under a weaker censorship assumption that Y and C are conditionally independent given X. Their method also demonstrated better small sample performance than the synthetic data EL method when the errors ϵ 's are independent and identically distributed (homogeneous). On the other hand, it is not easy to extend the method of Zhou and Li [18] to incorporate auxiliary information and to construct confidence regions for linear combinations of the regression coefficients, which are relatively easy using the synthetic data EL method (cf. Li and Wang [15]).

In this paper, we develop a new empirical likelihood method for linear regression with right censored data. We construct an estimated empirical likelihood based on the estimation equation for the Buckley-James [3] estimator. The approach inherits many appealing features of empirical likelihood. For example, the shape and orientation of a confidence region are entirely determined by data. Most importantly, it does not involve variance estimation. Our simulation shows that, under the assumptions of model (1), our method is far superior to the synthetic data empirical likelihood method of Li and Wang [15], with substantially shorter confidence intervals and higher coverage probabilities. Compared to Zhou and Li [18]'s method, our method may not be as efficient, but it is easier to compute. It can also be extended easily to incorporate auxiliary information and to construct confidence regions for linear combinations of the regression coefficients. More discussion of the pros and cons of our method in relation to existing methods are given later in Remark 2.

The use of empirical likelihood dates back at least to Thomas and Grunkemeier [19] who constructed confidence intervals for survival probabilities. The idea was later popularized by Owen [20, 21] who derived confidence regions for the mean of a random vector. Since the work of Owen [20, 21], there has been extensive developments of empirical likelihood methods for a wide range of applications including, among others, linear models [22-24], generalized linear models [25], quantile estimation [26], biased sample models [27], generalized estimating equations [28], dependent process model [29], partial linear models [30], mixture proportions [31], confidence bands for survival functions [32], confidence bands for quantile functions [33], censored cost regression [34], and confidence tubes for multiple quantile plots [35]. Some nice discussion of properties of empirical likelihood can be found in DiCiccio, Hall and Romano [36], Hall [37], and Hall and Scala [38], and others. A comprehensive survey of empirical likelihood and further references can be found in Owen [39] and Li et al. [40].

In Section 2, we derive an estimated empirical likelihood for β by combining the ideas of Owen [21, 22] and Buckley and James [3]. An adjustment factor is adopted so that the adjusted empirical likelihood has an asymptotic standard Chi-square distribution. We also discuss how to incorporate auxiliary information using empirical likelihood. Section 3 presents results from a simulation study to illustrate the performance of our method compared with the synthetic data empirical likelihood method. A real data example is also provided. Section 4 gives some concluding remarks. The proofs are collected in the appendix.

2. Empirical Likelihood Inference

2.1. Empirical Likelihood Based on the Buckley-James Estimating Equation

We motivate our procedure by first considering the case where μ = E(X) and F are known.

We first give a review of the Buckley-James equation. It can be shown that, under model (1), or, equivalently,

Because Y is not always observable, we impute Y by its conditional expectation given the observed data as follows:

Noting that E(Y * | X ) = E(Y | X ), we can replace Y by Y * in (3). This leads to where

Now, the problem of testing H 0 : β = β 0 is equivalent to testing based on n i.i.d. observations W (β 0, F, μ ), i = 1,…, n. This problem can be readily solved using the results of Owen [21]. Specifically, define Owen [21] showed that where λ is the solution of the equation Moreover, under H 0, l (β 0, F, μ ) has an asymptotic central Chi-square distribution with p degrees of freedom. Thus, one would reject H 0 if l (β 0, F, μ ) > χ 2, where χ 2 is the upper α quantile of the standard Chi-square distribution with p degrees of freedom.

Now, we consider the case where both F and μ are unknown. Define e = Z − X β, i = 1,…, n. We estimate F by the Kaplan-Meier [14] estimator of F based on {(e , δ ), i = 1,…, n}, where e (1) ≤ ⋯≤e ( are the order statistics of the e -sample, and δ ( is the δ associated with e ( , i = 1,…, n. In addition, we estimate μ by the sample mean .

We propose to use as a likelihood ratio statistic for testing H 0. However, we can no longer use Owen's [21] result for the null limiting distribution of because 's are not i.i.d. The following theorem states that an adjustment factor is needed so that the limiting null distribution is standard Chi-squared.

Theorem

Assume that the conditions listed in the appendix hold. Define that where

Then, under H 0 : β = β 0, where χ 2 is a standard Chi-square random variable with p degrees of freedom.

It follows immediately that an approximate α-level test rejects H 0 if Moreover, an approximate 1 − α confidence region for β is given by

Remark 2

Although both the above derived method and Zhou and Li [18] use the Buckley-James estimation equation, a sample version of (7), they are different in that we use the complete data likelihood, whereas Zhou and Li [18] uses the exact censored data likelihood to construct the EL. Similar to Li and Wang [15], the above method can be extended to incorporate auxiliary information and to obtain EL procedure for a subset, contrast, or linear combinations of the regression coefficients, which does not seem easy when using the method of Zhou and Li [18]. As an illustration, we show below how to extend our method to incorporate auxiliary information.

2.2. Empirical Likelihood with Auxiliary Information

Auxiliary population characteristics of the covariate X are sometimes available in practice. Effective usage of the auxiliary information can lead to more efficient inference (cf. Chen and Qin [41], Qin and Lawless [28] and Zhang [42, 43]). Here, we show how to use empirical likelihood to incorporate auxiliary information of X.

Assume that the available auxiliary information on X is given in the form Eg(X) = 0, where g(x) = (g 1(x),…, g (x)), r ≥ 1, is a vector of r known functions. To make use of the auxiliary information, we maximize subject to ∑ p = 1, ∑ p g(X ) = 0, and .

Let . By the method of Lagrange multipliers, it can be shown that (17) is maximized at where ζ satisfies the following equation Hence, the empirical log-likelihood ratio function for β is given by

Similar to the previous section, an adjustment factor is needed for l (β) to have a standard Chi-square asymptotic distribution, as stated in the following theorem.

Assume that V 1 = Eg(X )g (X ) is positive definite and that Eg(X)W (β, F, μ )) exists. Define that V (β) = (1/n)∑ g(X )g (X ), where and are defined in Theorem 1. Then, under the conditions of Theorem 1, as n → ∞, where and Ψ(β) = (∑ A (β))(∑ A (β)).

3. Numerical Results

We carried out Monte Carlo simulations to examine the performance of our proposed empirical likelihood method based on the Buckley-James estimating equation (ELEE) in comparison to the empirical likelihood method based on synthetic data (ELSD) [15, 16]. We considered five models. In model A, the data were generated from Y = 1 + X + ϵ, where X and ϵ are independent normal random variables with mean 0 and variances 0.25, respectively, the censoring time C is a normal random variable with mean μ and standard deviation 4. Model B is the same as model A except that X~ Bernoulli(0.5) − 0.5. Model C assumes that Y = X + ϵ, where X ~ N(0,0.52), ϵ ~ Weibull (shape = 1.843, scale = 1), and C ~ N(μ, 42). Model D is the same as model A except that C ~ N(μ + 2X, 15), allowing the censoring time to depend on X. Model E is the same as model A except that ϵ ~ N(0, X 2), allowing for heterogeneous errors. We adjust μ to produce different censoring rate (CR). We also vary the sample size n. The achieved confidence levels and average lengths of the ELEE and ELSD confidence intervals for the slope parameter are summarized in Table 1. Each entry in the table was computed using 3,000 Monte Carlo samples.

Table 1

Comparison of the coverage probability (CP) and average width (Width) of two empirical likelihood confidence intervals for the slope parameter under four different models with various sample size (n) and censoring rate (CR). Here, ELEE is the proposed method, and ELSD is the method of Li and Wang [15]. Each entry is based on 3,000 Monte Carol samples.

			Nominal level = 90%				Nominal level = 95%
Model	n	CR	CP		width		CP		Width
			ELEE	ELSD	ELEE	ELSD	ELEE	ELSD	ELEE	ELSD
	50	0.75	0.94	0.77	1.59	2.18	0.98	0.84	1.95	2.70
	100	0.75	0.94	0.82	0.85	1.66	0.97	0.89	1.03	2.02
	500	0.75	0.91	0.88	0.31	0.81	0.96	0.94	0.37	0.97
	50	0.3	0.94	0.87	0.69	1.30	0.97	0.92	0.82	1.55
A	100	0.3	0.93	0.89	0.45	0.94	0.97	0.94	0.53	1.12
	500	0.3	0.91	0.90	0.18	0.43	0.96	0.95	0.21	0.51
	50	0.1	0.95	0.88	0.59	1.10	0.98	0.93	0.70	1.30
	100	0.1	0.93	0.89	0.39	0.79	0.97	0.94	0.47	0.94
	500	0.1	0.90	0.90	0.16	0.36	0.95	0.95	0.19	0.43
	50	0.75	0.93	0.83	1.20	2.05	0.95	0.88	1.40	2.50
	100	0.75	0.94	0.87	0.77	1.49	0.97	0.92	0.92	1.80
	500	0.75	0.93	0.89	0.30	0.67	0.96	0.94	0.36	0.80
	50	0.3	0.95	0.88	0.66	1.23	0.98	0.94	0.78	1.48
B	100	0.3	0.94	0.90	0.44	0.88	0.97	0.95	0.52	1.05
	500	0.3	0.92	0.90	0.18	0.39	0.96	0.95	0.21	0.47
	50	0.1	0.94	0.89	0.58	1.08	0.97	0.95	0.69	1.29
	100	0.1	0.94	0.90	0.39	0.77	0.97	0.94	0.46	0.92
	500	0.1	0.91	0.91	0.16	0.35	0.96	0.95	0.19	0.41
	50	0.75	0.93	0.77	1.49	2.01	0.96	0.83	2.01	2.46
	100	0.75	0.93	0.82	0.80	1.56	0.97	0.88	0.97	1.90
	500	0.75	0.92	0.87	0.29	0.76	0.96	0.93	0.35	0.92
	50	0.3	0.93	0.86	0.67	1.21	0.97	0.92	0.81	1.45
C	100	0.3	0.93	0.88	0.44	0.88	0.97	0.93	0.53	1.05
	500	0.3	0.91	0.89	0.18	0.40	0.96	0.94	0.21	0.48
	50	0.1	0.94	0.87	0.60	1.01	0.97	0.93	0.71	1.21
	100	0.1	0.93	0.89	0.39	0.73	0.97	0.94	0.47	0.87
	500	0.1	0.92	0.90	0.16	0.33	0.96	0.95	0.19	0.39
	50	0.75	0.95	0.68	1.74	2.45	0.97	0.77	1.87	2.98
	100	0.75	0.94	0.60	0.83	1.83	0.97	0.69	1.01	2.21
	500	0.75	0.92	0.12	0.30	0.89	0.96	0.18	0.36	1.06
	50	0.3	0.94	0.81	0.68	1.31	0.97	0.88	0.82	1.56
D	100	0.3	0.93	0.76	0.45	0.94	0.97	0.84	0.53	1.12
	500	0.3	0.91	0.39	0.18	0.43	0.96	0.51	0.21	0.51
	50	0.1	0.94	0.86	0.59	1.09	0.97	0.92	0.71	1.30
	100	0.1	0.94	0.87	0.39	0.78	0.97	0.93	0.47	0.93
	500	0.1	0.91	0.78	0.16	0.36	0.95	0.86	0.19	0.42
	50	0.75	0.78	0.76	1.52	2.25	0.85	0.83	1.79	2.75
	100	0.75	0.78	0.80	0.85	1.74	0.85	0.87	0.62	2.11
	500	0.75	0.73	0.85	0.34	0.88	0.81	0.92	0.40	1.06
	50	0.3	0.81	0.85	0.76	1.43	0.87	0.91	0.89	1.71
E	100	0.3	0.79	0.87	0.53	1.05	0.86	0.92	0.62	1.26
	500	0.3	0.77	0.89	0.22	0.50	0.84	0.94	0.26	0.60
	50	0.1	0.81	0.86	0.69	1.24	0.87	0.92	0.81	1.48
	100	0.1	0.80	0.87	0.49	0.91	0.86	0.93	0.57	1.08
	500	0.1	0.76	0.89	0.20	0.42	0.83	0.94	0.23	0.50

Model A: Y = 1 + X + ϵ, where X ~ N(0,0.52), ϵ ~ N(0,0.52), and C ~ N(μ, 42); model B: Y = 1 + X + ϵ, where X∼ Bernoulli(0.5) − 0.5, ϵ ~ N(0,0.52), and C ~ N(μ, 42); model C: Y = X + ϵ, where X ~ N(0,0.52), ϵ∼ Weibull (shape = 1.843, scale = 1), and C ~ N(μ, 42); model D (Dependent censoring): Y = 1 + X + ϵ, where X ~ N(0,0.52), ϵ ~ N(0,0.52), and C ~ N(μ + 2X, 15); model E: Y = 1 + X + ϵ, where X ~ N(0,0.52), ϵ ~ N(0, X 2), and C ~ N(μ, 42).

We see from Table 1 that under models A–D, the coverage probabilities of ELEE method are consistently close to or slightly above the nominal level, whereas the ELSD method can have severe under-coverage for small samples (n = 50) and large censoring rate (75%). Furthermore, the ELEE confidence intervals are much narrower than the ELSD confidence intervals. In particular, the ELSD method failed completely with unreasonably low coverage under model D when the censoring time is dependent on X. On the other hand, under model E with heterogeneous errors and independent censoring time, ELEE showed larger coverage probability errors than ELSD, as one would have expected. Thus, the ELEE method seems to dominate the ELSD method when the errors are homogeneous, but can be outperformed by ELSD in the presence of heterogeneous errors.

We now illustrate our method using the Stanford heart transplant data (Miller [44], Table 1). The data include the lengths of survival (in days) after transplantation, ages at time of transplant, and T5 mismatch scores for 69 patients who received heart transplants at Stanford and were followed to April 1, 1974. Twenty-four patients were still alive on April 1, 1974 and thus their survival times were censored. For illustration purpose, we considered two models, labeled as (I) and (II), respectively, where the dependent variable Y is the logarithm to base 10 of the length of survival from transplantation. Specifically, model (I) regresses Y on the mismatch score T5, and model (II) regresses Y on age. As in Koul et al. [11], regression of survival on the mismatch score T5 was performed with nonrejection-related death being treated as censoring since the mismatch score is directed at the rejection phenomenon [44]. Table 2 reports the parameter estimates and 95% confidence intervals for the slope parameters using our empirical likelihood method based on the Buckley-James estimating equation (ELEE) and the empirical likelihood method based on synthetic data (ELSD) [15, 16].

Table 2

Empirical likelihood confidence interval estimates for heart transplant data. Nominal level = 95%.

Model		Parameter estimates		Confidence intervals
	Parameter	BJ	KSV	ELEE	ELSD
(I)	β T5	−0.593	0.258	(−2.740, 0.645)	(−0.596, 0.928)
(II)	β age	−0.028	0.055	(−0.065, 0.032)	(0.001, 0.128)

Note: BJ refers to the Buckley-James [3] estimate, and KSV refers to the synthetic data estimate of Koul et al. [11]. ELEE is the proposed empirical likelihood method based on the Buckley-James estimating equation, and ELSD is the empirical likelihood method based on synthetic data [15, 16].

It is seen from Table 2 that both the point estimates and confidence intervals are quite different between the ELEE and ELSD methods for this data. For example, the KSV slope estimates of T5 and age are positive, which seems to contradict to the common belief that the survival time tends to be negatively correlated with T5 (the mismatch score) and age at diagnosis. For model (II), the ELSD confidence interval does not include zero and thus concludes a significant age effect at the 5% significant level, whereas the ELEE method does not produce a significant result. The ELSD results could be misleading in this example since the censoring time appears to depend on age and T5 as pointed out by Leurgans [45] and Li and Lu [17]. We have also examined the Cox-Snell residual plots for a number of parametric accelerated failure time (AFT) models. We found that the log-normal AFT model fits the data fairly well which indicates that the semiparametric AFT model should also fit the data well.

4. Discussion

This research adds a new tool to the toolbox of empirical likelihood (EL) methods for linear regression with right censored data. The three EL methods, namely, the method of Li and Wang [15], the method of Zhou and Li [18], and the method of this paper should be regarded as complementary, as opposed to competing, methods for linear regression with right censored survival data. Each of the three approaches has its own merits and shortcomings. None dominates the others in every aspect. So it is important to understand the pros and cons of these methods. First of all, if the errors in the regression model are i.i.d (homogeneous) and the censoring time is conditionally independent of the survival time given the covariates, then the method of Li and Zhou [18] and the method of this paper are expected to be superior to the Li and Wang [15] method, as demonstrated by simulations in Table 1 and Li and Zhou [18]. This is because the former two methods use the Buckley-James estimating equations which take advantage of the homogeneous error assumption to implicitly impute a censored observation from the model using all residuals, whereas the latter uses synthetic data that does not utilize the homogeneous error assumption. Furthermore, based on our limited experience, the Li and Zhou [18] seems to be more efficient than the method of this paper. This is not a surprise since, Li and Zhou method [18] uses the exact censored likelihood, whereas the ELEE method of this paper uses an approximate likelihood. On the other hand, Li and Zhou [18] only developed an EL procedure for the whole vector of regression coefficients and did not discuss how to extend their method to incorporate auxiliary information which can be easily done using the approach of this paper as described in Section 2.2. Secondly, if the errors are heterogeneous, then the Li and Zhou [18] method and the method of this paper are expected to fail as indicated by Table 1 (model E), but the synthetic data method of Li and Wang [15] would still be asymptotically valid under the stronger censorship assumption that the censoring time is independent of both the survival time and the covariates. Finally, the method of Li and Wang [15] can fail completely when the censoring time is dependent on the covariate as shown in Table 1 (model D).

Some further extensions are warranted in future research. For example, the method of this paper can be extended to draw inference for linear combinations of the regression coefficients along the lines of Li and Wang ([15, Section 3]). It can also be extended to construct EL confidence regions based on estimating equations for a class of M-estimators described by Ritov [5]. Finally, this paper focuses only on EL methods for right censored data. It would be interesting to investigate how these EL methods compare to other methods such as the rank-based regression methods in future studies.