The survival data are common in many fields, such as economics, business, industrial engineering and biomedical studies. The nonparametric and semiparametric modeling has been studied extensively because it offers valid estimation and inference with less stringent model assumptions. The Cox proportional hazards model is the most used semiparametric regression models for the analysis of survival data (Cox, 1972) due to the availability of statistical software packages for implementation of the estimation and inference procedures. The accelerated failure time (AFT) models relates the logarithm of the failure time linearly to the covariates. The parameters in the model provides a direct interpretation. However, its estimation and inference procedures are challenging in the presence of censroing. Over the last forty years, there has been considerable research on the AFT models, see Buckely and James (1979), Jin et al. (2003, 2006), Koul et al. (1981), Lai and Ying (1991a, 1991b, 1995), Miller and Halpern (1982), Powell (1984), Prentice (1978), Ritov (1990), Robins and Tsiatis (1992), Tsiatis (1990), Wei et al. (1990), Yang (1997), Ying (1993), Zhou (2005a), Zeng and Lin (2007) and among many others. An earlier review on AFT model can be found in Wei (1992). In this paper, our goal is to review and discuss some newly developed practically useful statistical methods which specifically deal with the semiparametric accelerated failure time model.

Let T_i be the failure time for the i^th patient, i = 1, …, n. Due to censoring C_i, we observe Y_i = min(T_i, C_i) and δ_i = I{T_i ≤ C_i} which takes value 1 if T_i ≤ C_i and 0 otherwise. The semi-parametric AFT model is of form

where X_i, i = 1, …, n are observed covariates, β is p × 1 unknown parameter vector and ε_i, i = 1, …, n are independent and unobserved errors (the mean of the ε_i is not necessarily 0). The semi-parametric AFT model is essentially a regular semiparametric linear regression model and its regression coefficient β is directly quantifies the impact of covariates X on the survival time instead of the more abstract hazard rates. Note that the logarithm in the model may be replaced by any other known strictly increasing functions. However, in the presence of censoring, the intercept parameter is not included in the model specification because the intercept parameter cannot be estimated well as pointed out in Wei (1992).

It is assumed that conditional on the p × 1 covariate vector X_i for the i^th subject, the censoring variable C_i and the failure times T_i, i = 1, …, n are independent, and ε_i, i = 1, …, n are independent and identically distributed random variables whose common distribution function F(·) is completely unspecified. (Cox and Oakes, 1984; Kalbfleisch and Prentice, 2002; Lawless, 2003).

Let ei(β0)=log Yi-XiTβ0=min{log Ti-XiTβ0,log Ci-XiTβ0}=min{ɛi,log Ci-XiTβ0}, and N_i(β; t) = Δ_iI{e_i(β) ≤ t}, where I(·) is the indicator function that takes value 1 when the condition is satisfied 0 otherwise.

3.1. Gehan-type rank estimator

Currently, all available estimation method starts with following estimation equation, based on the Gehan’s statistic:

The function UnG is component-wise monotone in β (Fygenson and Ritov, 1994). It was shown that solving U_n(β) = 0 or ||U_n(β)|| = 0, is equivalent to minimizing G_n(β) (Jin et al., 2003), where

where the notation a⁻ means the negative part of a (e.g. 3⁻ = 0, (−2)⁻ = 2). The minimization of G_n(β) can be done by following linear programming approach: Minimize ∑i=1n∑j=1nΔiuij under linear constraints

Or alternatively, the minimization of G_n(β) is equivalent to minimization of

where M is the prespecified extremely large number (Jin et al., 2003). The minimization can be done with many available statistical softwares, such as the function l1fit in R.

Therefore, point estimator of β can be obtained, which is denoted as β̂_G. Under regularity conditions, it was shown that β̂_G is consistent and n(β^G-β)~N(0,Γ) as n → ∞, where Γ=AG-1BGAG-1, the A_G and B_G are defined in the expression (3.6) and (3.7) with φ(β,t)=∑j=1nI{ej(β)≥t}.

3.2. Weighted logrank-type rank estimator

The general weighted log-rank estimating function for β takes the form

or

where φ is a possibly data-dependent weight function. When φ = 1, the resulting U_φ(β) corresponds to the log-rank (Mantel, 1966) statistics. When φ(β,t)=∑j=1nI{ej(β)≥t}, the resulting U_φ(β) correspond to the UnG(β) in (3.1), which is Gehan (1965) statistics.

Under regularity conditions, the root β̂_φ of the estimating function U_φ(β) satisfies that n^1/2(β̂_φ − β) is asymptotically zero-mean normal with covariance matrix Aφ-1BφAφ-1, where

A^⊗2 = AA^T for a matrix A (Ying et al., 1992), λ(·) is the common hazard function of the error terms, and λ̇(t) = dλ(t)/dt (Lai and Ying, 1991b; Tsiatis, 1990; Wei et al., 1990; Wei, 1992; Ying, 1993).

In general, it is difficult to solve the equation U_φ(β) = 0 because U_φ(β) is neither continuous nor componentwise monotone in β. Jin et al. (2003) overcame the difficulty by a class of monotone estimation functions that approximate the general weighted log-rank estimating functions (3.4). The idea is to use the Gehan-type rank estimator β̂_G in the Section 3.1 and an iterative algorithm. Specifically, Jin et al. (2003) proposed to use following modified estimating function of (3.4):

where ψ(b; x) = φ(b; x)/S⁽⁰⁾(b; x), S(0)(β;t)=∑i=1nI{ei(β)≥t} and β̂ is a intial consistent estimator of β₀, say β̂_G. It is easy to see that

The Ũ_φ(β;β̂) is monotone in each component of β and is the gradient of the following function:

Similar to (3.2), the minimisation of L_φ(· ;β̂) can be implemented via linear programming.

The iterative algorithm for estimating β starts from setting: β̂₍₀₎ = β̂_G, and β̂_(k) = arg min_βL_φ(β; β̂_(k−1)) (k ≥ 1). If β̂_(k) converges to a limit as the number of iterations k → ∞, then the limit must satisfy U_φ(β) = 0. Therefore, the root β̂_φ of original estimating equation U_φ(β) = 0 can be obtained by the iterative algorithm. Under regularity conditions, Jin et al. (2003) showed that each estimator obtained in the k iteration β̂_k is a legitimate estimator of β (consistent and asymptotically normal) when the initial estimator is consistent and asymptotically normal, and if the β̂_k converges as k → ∞, the limit is indeed the solution of the U_φ(β) = 0. Specifically, Jin et al. (2003) showed that β̂_k is asymptotically a weighted average of β̂_φ and β̂_G,

where I is the p × p identity matrix,

and ψ̇₀(t) is the derivative of ψ₀(t), the limit of ψ(β₀; t) as n → ∞.

3.3. Least-squares-type estimator

When there is no censoring, T_i, i = 1, …, n are completely observed, the classical least-squares estimator of β is obtained by minimizing ∑i=1n(log Ti-XiTβ)2 in terms of β. The minimizer is the solution of the equation

where X¯=n-1∑i=1nXi.

In the presence of censoring, we can only observe Y_i and Equation (3.13) cannot be used directly. Buckly and James (1979) proposed to replace log T_i by following

where F̂_β(t) is the Kaplan-Meier estimator of F(t) based on {e_i(β), δ_i} (i = 1, …, n), i.e.,

The resulting Buckly-James estimator is the solution to the follwoing equation:

The difficulty to solve the Equation (3.15) lies in that the function ∑i=1n(Xi-X¯)(Y^i(β)-XiTβ) is neither continuous nor componentwise monotone in β. To overcome these difficulties, Jin et al. (2006) developed an iterative procedure to obtain a class of consistent and asymptotically normal estimators. Let Y¯(b)=n-1∑i=1nY^i(b), and define

or equivalently

Set U(β, b) = 0, which yields

It leads to an iterative algorithm

With the consistent initial estimator β̂₍₀₎ = β̂_G, the Gehan-type estimator in the Section 3.1, Jin et al. (2006) showed that β̂_m is consistent and asymptotically normal for every m ≥ 1, and β̂_(m) is asymptotically a weighted average of the Buckley-James estimator β̂_BJ and β̂_G,

where D=limn→∞(1/n)∑i=1n(Xi-X¯)⊗2 and A is the slope matrix of the Buckley-James estimating function.

3.4. Kernel-smoothed profile likelihood estimator

Zeng and Lin (2007) proposed a kernel-smoothed profile likelihood estimator by an approximate nonparametric maximum likelihood method. Based on the semiparametric AFT model (2.1), the observed data (Y_i, δ_i, X_i), i = 1, …, n yields the following log-likelihood function of β and hazard function λ(·) of the error ε:

Zeng and Lin (2007) showed that the direct maximization of the log-likelihood does not yield a solution. It also does not lead to the maximization of a profile likelihood obtained by replacing the λ(e_i(β)) by an estimator δi/∑j=1nI{ej(β)≥ei(β)} nor the maximization of a sieve profile likelihood obtained by replacing the λ(·) by a piecewise constant approximation. On the other hand, an approximation of

by

where K(·) is a kernel function (such as the standard normal density function) and a_n is a bandwidth, yields a kernel-smoothed approximation of the likelihood of β,

The kernel-smoothed profile likelihood estimator β̂_pl is obtained by maximizing the lns(β) and its variance-covariance estimator is obtained by the use of the information matrix of the lns(β).

Given β̂_pl, the λ(·) and Λ(·) can be estimated by

Under regularity conditions, Zeng and Lin (2007) showed that the estimator β̂_pl achieves semiparametric efficiency. The challenge for the practical use is the choice of the Kernel function, bandwidth, and the stability of the lns(β). Currently, there is no statistical software that implements the procedure.

Except the kernel-smoothed profile likelihood estimator, the variance-covariance estimation of the Gehan-type rank estimator, weighted logrank-type estimator, and least-squares-type estimator is difficult because the corresponding variance-covariance matrices involve nonparametric estimation of the underlying probability density function. Here two computationally feasible approaches will be presented: resampling method and induced smoothing method.

4.1. Resampling method

The resampling method is similar to the resampling scheme similar to those in Rao and Zhao (1992), Parzen et al. (1994) and Jin et al. (2001). The basic idea is to perturb the objective functions that are minimized or estimating equations that are solved. Let Z_i (i = 1, …, n) be n independent positive random variables satisfying E(Z_i) = Var(Z_i) = 1 and be independent of the observed data (T̃_i, δ_i, X_i) (i = 1, …, n). The Z_i can be generated by a number of known distributions, such as the standard exponential distribution or Poisson distribution Poisson(1).

4.2. Induced smoothing method

Brown and Wang (2005) proposed a general variance estimation procedure based on an induced smoothing for non-smooth estimating functions. The approach is computationally efficient and easy to implement when the smoothing in terms of integration with respect to a Gaussian distribution has an explicit form. For the Gehan-type rank estimator, the induced smoothing in Brown and Wang (2005) yields an explicit form and its varaince-covariance matrix can be easily estimated as shown in Brown and Wang (2007). However, the induced smoothing in Brown and Wang (2005) for the general logrank-weighted estimator and least-squares-type estimator does not yield an explicit form and cannot be used directly. To overcome the difficulty, Jin et al. (2015) developed a general Monte carlo approach for variance estimation for estimators based on induced smoothing which can be applied to both the general logrank-weighted estimator and least-squares-type estimator.

In this paper, we have reviewed recently developed methods for point and variance estimation for the AFT model in the analysis of right censored data. The development is based on the assumption that errors are independent and identically distributed. The resampling approach offers valid inference but numerically intensive for large datasets, on the other hand, the induced smoothing approach is attractive due to its computational efficiency.

The paper of Jin and Ying (2004) studied asymptotic theory of rank estimation for AFT model under fixed censorship. Zhou (1992) and Jin (2007) studied M-estimation for the AFT models.

For the hypothesis testing framework, empirical likelihood approach has also been developed, see Zhou (2005a, 2005b), Zhou and Li (2008). For the heteroscedastic errors, Stute (1993, 1996) studied convergence properties of weighted estimators and Zhou et al. (2012) developed an empirical likelihood approach under a hypothesis testing framework.

For time-dependent covariates, there has been theoretical development, see Robins and Tsiatis (1992) and Lin and Ying (1995). The profile nonparametric likelihood approach of Zeng and Lin (2007) is applicable theoretically, but there has been no numerical investigation.

There are still many research problems, such as model checking, variable selection and efficient and reliable software development.