Buckley-James Estimator of AFT Models with Auxiliary Covariates

In this paper we study the Buckley-James estimator of accelerated failure time models with auxiliary covariates. Instead of postulating distributional assumptions on the auxiliary covariates, we use a local polynomial approximation method to accommodate them into the Buckley-James estimating equations. The regression parameters are obtained iteratively by minimizing a consecutive distance of the estimates. Asymptotic properties of the proposed estimator are investigated. Simulation studies show that the efficiency gain of using auxiliary information is remarkable when compared to just using the validation sample. The method is applied to the PBC data from the Mayo Clinic trial in primary biliary cirrhosis as an illustration.


Introduction
It is not uncommon to have one or more missing or mismeasured covariates in large cohort epidemiological studies. There are always cases in medical studies, where it is difficult to obtain an accurate measurement for all patients due to a procedure being too expensive or invasive. Alternatively, some auxiliary measurements which are less precise, but highly related to the target procedure, can be easily collected. In some situations, all of the measurements are error prone, while in other cases, a validation subsample, where the measurements are all accurately taken, is made available.
The former is a pure measurement error problem. The purpose of this paper is to investigate the inference of the latter cases in a failure time setting. In some cases, the validation sample could be large enough on its own, so one could choose to ignore all data from subjects that have missing or mismeasured values for any of the covariates, with just a minor efficiency loss. However, if the validation sample is relatively small, utilizing the auxiliary information will lead to remarkable efficiency gain, as our simulation results will show.
However, due to the direct physical interpretation of the AFT models, and the fact that AFT models are robust to model misspecification in the sense that ignoring a covariate will not lead to much bias in estimating the remaining regression coefficients, see [7] (Cox and Oakes, 1984), the biasing effect of covariate measurement error on AFT models deserves further investigation. A recent work on the subject of measurement error in AFT models was done by [16] (He et al., 2007), using a simulation and extrapolation approach. [17] (Yu and Nan 2010) studied the regression calibration approach within the semiparametric framework, assuming a known parametric relationship between the accurately measured covariates and their auxiliaries, up to a few unknown nuisance parameters. [18] (Granville and Fan, 2012) studied the parametric AFT models with auxiliary covariates based on maximum likelihood method.
In this paper, we study the Buckley-James estimator [19] (Buckley and James, 1979) of AFT models with auxiliary covariates. The Buckley-James estimator was shown by [20] (Jin et al., 2006) to be consistent and asymptotically normal when using a consistent estimator as the initial value, due to its asymptotic linearity. Some other insights about the consistency and asymp-totic theory of this estimator has been investigated by [21] (James and Smith, 1984) and [22] (Lai and Ying,1991), among others. We propose a local polynomial approximation method to handle the missing or mismeasured covariates, through the estimation of the conditional expectation of the unobservable estimating functions. This approach makes neither distributional assumptions about the model error term e i , beyond it having mean zero and a finite variance, nor parametric assumptions on the relationship between the correctly measured covariates and their auxiliary variables. The proposed approach will be introduced through a kernel smoothing method, a special case of the local polynomial approximation, see [14] (Fan and Wang, 2009), mainly due to the ease of presentation. See [23] (Nadaraya, 1964), [24] (Watson, 1964), and [25] (Wand and Jones, 1995) for details of kernel smoothing. Intensive simulation studies were conducted to investigate the small sample performance of our proposed method. The results show a remarkable efficiency gain over the method which ignores the auxiliary information.
The remainder of this paper is organized as follows. In the second Section, we introduce Buckley-James estimator for the accelerated failure time model and present our estimation method. Then we investigate the asymptotic properties of our proposed estimator. The Section thereafter contains the results and discussion of our numerical studies, including simulations and the PBC data illustration. In the last Section, we put forth some concluding remarks. The proofs for Theorems were deferred to the appendix.

Inference Methods of Accelerated Failure Time Model
Let T i and C i , i~1, . . . ,n be the failure and censoring times for the ith subject in a large study cohort. Due to the censoring, we observe S i~m in (T i ,C i ) as well as a failure indicator d i~I (T i ƒC i ). Let fX i ,Z i g denote the covariate vector where X i is the component which is only observed in the validation set and Z i is the component that is available for the full study cohort. Let W i be the auxiliary covariate to X i . Hence the data consists of the validation sample fS i , d i , Z i , X i , W i g, and the nonvalidation sample fS i , d i , Z i , W i g. In this paper we assume that X i is scalar, mainly because of the simplicity of the presentation, and Z i may be either a scalar or a vector. In practice, X i could also be closely correlated with Z i . A special case is the classical measurement error model W i~Xi zU i , where U i is the error encountered when measuring X i . It is assumed that the U i 's are independent and identically distributed random normal variables, U i *N(0,s 2 u ). Of the n observations, the validation sample contains n V observations, and the non-validation sample contains n V V~n {n V observations.
The accelerated failure time model based solely on the validation sample, can be expressed as where b'~(b 1 ,b' 2 ) is a vector of unknown regression coefficients and the e i 's are independent and identically distributed with an unspecified distribution F which has mean zero and finite variance. Equation (1) assumes automatically that W i provides no additional information about the failure time, given fX i , Z i g.
Without making any assumption to the distribution of e i , the Buckley-James least squares procedure (Buckley and James, 1979) estimates the regression parameters through the minimization of The least squares estimates of b 1 and b' 2 are such that and In Then and However, the distribution of e i is unknown. The distribution of Y i , and consequently, E½Y i DY i w log (C i ) are both unknown. The censored observations are hence replaced bỹ . . ,n are the residuals, and is the Kaplan-Meier Product Limit estimator of the distribution function of the residuals, F . c F b F b (e) is a discrete function which will not tend to 1 as e increases if the largest residual is censored. Therefore, following the convention of Buckley and James, the largest residual is redefined as uncensored for all calculations, if necessary. LetỸ is the solution of the following equation, It should be noted that c(b) depends on X i , which is available only for the validation sample. For the non-validation sample we can substitute the estimates of their conditional expectations given the auxiliary and other available covariates. The local polynomial approximation approach can be applied for this purpose, see [14] ( Fan and Wang, 2009). For the simplicity of the presentation, we use the kernel smoothing method to estimate the conditional expectation of the unobserved covariates given the auxiliary information.
Note that this simplification does not necessarily lead to efficiency loss. Since the direct estimation of the conditional expectation of the estimating function depends also on the Kaplan-Meier estimation of the survival function of the regression residuals, it could also introduce additional instability into the inference, as compared with imputing the estimated conditional expectation of the mismeasured covariate. Our simulation also revealed this observation (results not included).
The conditional expectation of the mis-measured covariate, where is the chosen vector of bandwidth. Using these b X i X i 's in place of the n V missing X i 's, we may solve (6) for b using a numerical method, like Broyden's method. Note that, Broyden's method requires two initial values, while the method of [21] (James and Smith, 1984) only requires a single initial value. In this function, c n (b (k) ) is the value of (6) when calculated using b (k) . A very natural selection of the initial value of b is the least squares regression estimator calculated from the validation sample.
The standard deviation of these estimators is estimated using bootstrapping. For R replicates, a simple random sample with replacement of the full sample size is taken from the observed data and the above estimation method for b is repeated on each replicate. A sample standard deviation is then calculated to estimate the true standard deviation of the b estimators.
Remark 1 In order to retain the same quality of information among the replicated estimations, an alternative method of resampling was attempted to keep the proportion of censored observations constant in each replication. We defined d Ã~Pn i~1 d i to be the total number of observations with uncensored failure times. For R replicates, a simple random sample with replacement of size d Ã was taken from these uncensored observations, and a simple random sample with replacement of size n{d Ã was taken from the remaining censored observations. However this alternative method was found to underestimate the true standard deviations and resulted in coverage probabilities that were lower than the nominal level. The reason of this outcome is mostly due to the fact that the independence of the censoring mechanism was broken by the sampling method.

Defining a Solution
In order to solve the estimating equations for the regression parameters, we use the iterative scheme of b (kz1)~cn (b (k) ). However, as noted by [19] (Buckley and James, 1979) and [21] (James and Smith, 1984), these iterations need not converge. The c(b) function is discontinuous and piecewise linear in b so an exact solution may not exist. When this is the case, the iterations can oscillate between two values of b. We define a possible alternate solution which is closest to satisfying b~c(b), or c(b){b~0. If b (k) is oscillating between two points due to the lack of an exact solution, we define the alternate solution as b (k) that minimizes the modulus of this difference, When the iterations do not converge, a cut-off point has to be determined to stop the iterations. It is advised to select a number of iterations that is slightly greater than the amount required for the convergence when a solution typically does exist. For most of our simulations, we have set this point at k~11. In many cases, our simulations converge in three or four steps. So at k~5 or k~6, b (k) breaks the loop when checked against b (k{1) for convergence, implying the solution being reached at iteration k{1. If the iteration does not converge, the first ten values are checked and whichever value, after, say 5 steps of iteration, satisfies equation (8) is selected as a solution.
When dealing with real data, it is advised to choose an arbitrarily large number for the cut-off point to find the best possible solution.

Asymptotics
In this section, we investigate the asymptotic properties of our proposed estimator. For that sake, we rewrite the estimating function and the Kaplan-Meier estimator of the residual survival function in the counting process frame work. Define a function The estimating equation (6) can be rewritten as The Buckley-James estimator solves the above equation. When some of the covariates are accurately recorded only for the validation sample, but with relevant auxiliary information available for the whole study cohort, the estimating functions involved those mis-measured covariates belonging to the nonvalidation sample. We propose to estimate those terms by using the local polynomial approximation approach. Let n denote the size of whole study cohort, g i , for i~1, . . . ,n be the validation The derived estimating equation is then Our proposed estimator of the regression parameter, accommodating the auxiliary information, b b b solves this derived estimating equation.
For a vector a, define a 60~1 , a 61~a , a 62~a a' and Without loss of generality, let d be the dimension of C i in the definition of the local polynomial approximation. Suppose further that a is the order of the kernel function K, i.e.
The bandwidth conditions are given below.
The following assumptions, beyond the bandwidth conditions, are necessary for the asymptotic properties of the proposed method. C.0 The hazard rate function l(t) of e i (b) is such that Ð ? {? l(t)dtv?. C.1 There exists a constant Bw0, such that EZ i EƒB, EX i EƒB and EW i EƒB for all i~1, . . . ,n. C.2 F has a twice-continuously differentiable density f such that

Remark 2
The assumption C.3 is proposed just to simplify the proof of the asymptotics of the Buckley-James estimator. This could be violated due to the instability of the Kaplan-Meier estimator of the survival function when getting into the distribution tail. When this violation happens, the tail modification by [22] (Lai and Ying, 1991) should be applied and the method of [26] (Jin, et al., 2006) of selecting a consistent and asymptotically normal estimator as the initial value can be adopted. U U(b) converges in distribution to a zeromean normal random vector with covariance matrix rS(b)z(1{r)S 1 (b), where r~lim n?? n V =n, where Q~lim n??
and S 1 is defined as Theorem 2 Under assumptions C.0-C.4 and the bandwidth conditions is asymptotically linear within the n {1=3 neighborhood of b, with probability 1, in the sense that and A Ã as

Proof of Theorem 1
Let V be the set of the indices of the validation sample, V V that of the non-validation sample and g i~I (i[V ), i~1, . . . ,n be the validation indicator. Let Let Further, and  By the martingale representation of Kaplan-Meier estimator of the survival function see [25] (Fleming and Harrington, 1991), also see [20] (Jin, et al., 2006), the bandwidth condition [BC] and the continuity of F (assumption C.2), the first two terms in the above equation can be rewritten as and The Kaplan-Meier estimators of the survival functions of b e e i (b) and e i (b) lead to and ð ? and This term can be further rewritten as By assumptions C.0, C.2 and Lenglart's inequality, we have Hence the estimating function can be rewritten as is a sum of n V i.i.d. terms hence central limit theorem applies. By conditions C.0 through C.4 and the martingale central limit theorem, U Ã V (b) converges in distribution to a normal random vector. Further, by independence ofŨ

Proof of Theorem 2
From the equation (U) in the proof of Theorem 3.2, and by Theorem 4.1 of Lai and Ying (1991), we have, for DDb{bDDvn {1=3 , with probability 1. The termŨ U V (b) consists of two parts, we have with probability 1. Hence with probability 1. Corollary 1 and Theorem 3 are direct conclusions of Theorems 1 and 2.

Simulation Studies
In this section we examine the small sample performance of our proposed estimator. Let b b b S denote our proposed estimator of the regression coefficients. Its small sample performance is compared with three alternative estimators: the validation estimator ( b b b V ) which is based solely on the validation sample; the naive estimator ( b b b N ), which ignores the measurement error by assuming that the unobserved X i 's are equal to the observed W i 's; and the complete case estimator ( b b b CV ), when we assume that X i are observed for the whole study cohort.   Table 2. Results after 500 simulations for b'~( log (2), log (1:5))~(0:693,0:405) using an extreme value error term.
n n V Censor Rate The data for the simulations was generated in the following way. The X i and Z i are generated from a uniform distribution, X i ,Z i *uniform½0,5. For each X i , the auxiliary covariate is defined as W i~Xi zU i , where U i is generated from a normal distribution with mean zero and standard deviation s u . The value of s u determines the magnitude of the measurement error. The failure times were then defined as T i~e xpfY i g where Y i~b1 X i zb' 2 Z i ze i . The e i 's were taken to be independent and identically distributed from either a standard normal, standard extreme value, or logistic distribution, respectively.
Various other parameters are controlled over all simulations. Each run calculates 1000 replicates in the bootstrapping to give consistent estimators of the standard deviations. The parameters were chosen as b'~(b 1 ,b 2 )~( log (2), log (1:5)). Within a simulation, the censoring times are randomly generated from a uniform distribution with lower limit 0 and an appropriate upper limit to ensure an approximate 30% or 50% censor rate. The n and n V values are chosen to be either n~400 and n V~2 00, having half of the data in the validation set, or n~250 and n V~1 50, with the validation set containing 60% of the data. Finally, two values of s u are selected, s u~0 :5, and s u~0 :8. For the kernel smoothing used to calculate b b b S the Gaussian kernel function is selected, which has an order of 2, where u~(W i {W j )=h. We choose bandwidth h~2s u n {1=3 as used by [12] (Zhou and Wang, 2000).
The standard error (SE), standard deviation (SD), and coverage probability (CP) are calculated for each set of simulations. The SE values are the sample standard deviations of the b estimates, the SD values are the mean standard deviations generated from the bootstrapping in each simulation, and CP is equal to the percentage of simulations that had the true b value within a 95% confidence interval around its estimate when using the result of the bootstrapping for the standard deviation. The results are presented in Tables 1, 2, and 3.
From Tables 1, 2, and 3, we make the following observations: The coverage probabilities for the 95% confidence intervals are very close to their nominal level, except for b b b N when s u is large, where the estimate is severely biased.
(vii) The model experiences the least variation when the error term follows the standard normal distribution, with Table 4. AFT model analysis of PBC data, smoothing for log (ast).  Table 5. AFT model analysis of PBC data, smoothing for log (copper).

Application to PBC Data
To illustrate how to use the smoothing method in practice, we analyze the data from the Mayo Clinic trial in primary biliary cirrhosis (PBC) of the liver. PBC is a chronic liver disease that inflames and slowly destroys the bile ducts in the liver, impairing its ability to function properly. It is believed to be a type of autoimmune disorder where the immune system attacks the bile ducts. PBC occurs primarily in women, with approximately 90% of patients being women, most often between the ages of 40 and 60. There is currently no known cure for the disease; the only known way to remove PBC is through a liver transplant, see [27,28].
In the Mayo Clinic trial, 418 patients were eligible. Of these patients, mostly complete data was obtained from the first 312 patients. The remaining 106 patients were not part of the actual clinical trial but had some basic measurements taken and were followed for survival. The variables we used in our regression on the logarithm of time were age, patient's age (in years); albumin, serum albumin (in mg/dl); ast, aspartate aminotransferase (in U/ ml), once referred to as SGOT; bili, serum bilirunbin (in mg/dl); copper, urine copper (ug/day); edema, equal to 0 if no edema, 0.5 if untreated or successfully treated, or 1 if there exists edema despite diuretic therapy; protime, standardized blood clotting time. Of these, two cases were examined using either ast or copper for our X covariate to be smoothed due to incomplete data, while the others are mostly complete and thus are included in Z.
Edema was split into two categorical variables, edema05 and edema1, defined as & We also took the log transformation of albumin, ast, bili, copper, and protime, in the interest of making their marginal distributions closer to normal. For the smoothing of the unobserved log (ast) and log (copper) values, log (bili) was chosen as the auxiliary covariate for both due to its high correlation (w0:5) with both variables. The bandwidth h was calculated using the sample standard deviation of log (bili), resulting in s u &1:020551 and h~0:2734221 using the same formula as in the numerical simulations, h~2s u n {1=3 . Any observations missing a value for either of the Z covariates were removed, leaving both cases with n~416 while n V~3 12 for the model using log (ast) and n V~3 10 for the model using log (copper).
Examining Tables 4 and 5, we see that both X variables had their estimated standard deviations increase by a small amount due to the error added into the model from smoothing for a missing covariate instead of a mismeasured one. If W was of the form W i~Xi zU i like in the simulations, it could have resulted in a higher correlation between the auxiliary variable, W , and the X variable, depending on the magnitude of the measurement error. Despite the small increase in standard deviation the log (ast) term becomes significant at the 5% confidence level after smoothing, and while log (copper) was already significant, the p-value did decrease. For the Z variables, we see that they all have a smaller or approximately equal standard deviation after smoothing, which is expected when using the full sample size without needing smoothing for those variables, except for log (protime) and the intercept term which increased.

Discussion
In this paper we proposed the use of the Buckley-James estimator as a nonparametric method of estimating the regression parameters of an accelerated failure time model with auxiliary covariates. Kernel smoothing was applied using the auxiliary covariates to estimate missing or mismeasured covariates. The Buckley-James method is then applied to the whole study cohort for the inference of the covariates effect. The standard deviations of the estimates of the regression coefficients are estimated through bootstrapping. The proposed estimator is consistent and asymptotically normal.
This method was most effective in the case of mis-measured data due to the naturally high correlation between the corresponding W and X ,Z variables which resulted in the estimator involving the smoothing, b b b S , being more efficient than the validation estimator, b b b V , as shown in the numerical simulations. The method should also perform well for the missing variable case given a sufficiently strong correlation. The method was applied to the PBC data as an illustration.
The smoothing model is set up in a general format. In applications, we should only choose those variables which are highly related to the mismeasured one. By doing so we can avoid the situations such as the auxiliary covariates only occupy a narrow region, which could cause instability in the local smoothing, hence the whole model. Caution should also be taken when the proposed method is applied to a data with extremely small validation sample. A classic measurement error model might be a better option, where one can estimate the measurement error variance using the validation sample.