Weighted estimation of the accelerated failure time model in the presence of dependent censoring.

Independent censoring is a crucial assumption in survival analysis. However, this is impractical in many medical studies, where the presence of dependent censoring leads to difficulty in analyzing covariate effects on disease outcomes. The semicompeting risks framework offers one approach to handling dependent censoring. There are two representative estimators based on an artificial censoring technique in this data structure. However, neither of these estimators is better than another with respect to efficiency (standard error). In this paper, we propose a new weighted estimator for the accelerated failure time (AFT) model under dependent censoring. One of the advantages in our approach is that these weights are optimal among all the linear combinations of the previously mentioned two estimators. To calculate these weights, a novel resampling-based scheme is employed. Attendant asymptotic statistical results for the estimator are established. In addition, simulation studies, as well as an application to real data, show the gains in efficiency for our estimator.

Introduction

In medical studies, it is very common that death or withdrawal of study and progression on disease of interest simultaneously occur in the study. For this case, death or withdrawal of study may censor the development of disease. This type of data structure is called ‘semicompeting risks data’ [4].

Semicompeting risks data have been widely studied in the past decade. Some researchers used a Gamma copula to estimate the association parameter between the event of interest and dependent censoring [2], [4]. There is a literature that extended the methodology of [2] to the case that a nuisance parameter exists and also considered a more general copula model [22].

On the other hand, other researchers used semiparametric regression to model the event of interest and dependent censoring jointly. One approach is an estimation procedure based on the accelerated failure time (AFT) model [13], [17]. They used the artificial censoring technique to adjust the bias of the usual estimator. While the estimating equation of [13] is a U-statistic of order one, that of [17] is a U-statistic of order 2.

However, none of these papers fully discussed optimality of the estimator. In this case, choosing an estimator that is optimal from an efficiency viewpoint is an important issue for consideration. Here, we adapt the idea of [25], which proposed an optimal estimator whose form is a linear combination of estimators for multivariate failure time data. They used idea of [24], which proposed using combinations of dependent tests in the presence of missing values. Idea of [24] is to create a test which can maximize power based on linear combination of test statistics. Approach of [25] is simple and flexible, so it is sensible to apply their method in our case.

In this paper, we propose a weighted estimator by using methodology from [25]. Our weighted estimator combines those of [13] and [17]. The structure of this paper is as follows. In methods section, we review estimators proposed by [13] and [17] briefly. In addition, we describe details on our new weighted estimator. In model checking section, model checking procedure is briefly discussed. In simulation studies section, results of simulation studies will be given. Application of our method to a real data example is presented in real data analysis section. Some discussion concludes discussion section.

Methods

Review of Model

Let X be time to the event of interest, D the time to dependent censoring and C the time to independent censoring. All these times are transformed on a logarithmic scale. Let and . Define , Δ = I(D ≤ C) and let Z be covariates. The data contain n independent and identically distributed observations . The model is where 0 and 0 are k × 1 vectors, and are error terms with an unknown joint distribution. In this case, we assume that the model is identifiable only in upper wedge X < D [4], [17]. We assume that ϵ has unknown distribution H. The goal is to obtain an unbiased estimator of = ( , ) without nonparametrically estimating the distribution of ϵ , i = 1,…, n. We further assume that given Z, C and (X, D) are independent, but X and D can be dependent given Z. Now we are going to describe the procedures of [13] and [17] in turn.

Since only depends on independent censoring, a standard rank regression approach is available for estimation [11], [13], [15], [17], [20], [23]. The estimating equation for is given by where . The estimator of can be obtained by solving S () = 0.

For estimation of , simply replacing to does not yield unbiased estimation of . This is because the cause-specific hazard function for depends on , which violates the model assumption [13]. To fix this problem, many authors use artificial censoring techniques [3], [6], [7], [10], [13], [17]. In [13], a single constant term g() is proposed so that the estimation equation will be unbiased for estimation of in the two-sample case. The form of g() is . The proposed estimator in [13] is obtained by solving = 0, where and a∧b means minimum of a and b. In [17], pairwise comparisons of all the subjects is proposed so that each subject has different degree of the artificial censoring. The transformations suggested by [17] are The proposed estimator according to [17] is obtained by solving , which is defined by Note that X and D are not observable, but we can express transformation (Eq 1) and (Eq 2) by using observable quantities [7], [13].

Weighted estimator

Given these two estimation procedures, it is natural to consider their efficiencies with respect to standard error. However, in this point of view, neither estimator is superior to the other. Moreover, these estimators may not be optimal estimators with respect to the standard error. There is an argument that estimator of [17] gains more efficiency than that of [13] because pairwise comparisons lead to less artificial censoring than that in [13]. However, this logic only holds when we look at performance of estimators in the view of bias and variance across the estimators in simulation study. Concentrating on standard error of an estimator in a single dataset, the estimator by [17] may not provide better estimator than that of [13]. This will be seen in the real data analysis section.

The reason for this is due to estimation procedure of [17]. As discussed [7], for n samples, the number of comparisons of [13] for artificial censoring is of order n, while that of [17] is of order n 2. By definition of g (), different degrees of artificial censoring is applied to observations. It may lead more variation between observations, which makes standard error larger than that of [13].

Having discussed our data structure and estimators from [13] and [17], we now describe the proposed estimation in this paper. Let be estimator of 0, be estimator of 0 by [13] and be estimator of 0 by [17]. and are asymptotically unbiased estimators of 0.

We extend the scope of estimators which provide consistent estimation of 0. The natural extension of estimators of [13] and [17] is to consider collections of estimators that are linear combination of these two estimators with sum of weights being 1. By choosing proper weights, we can expect that the variance of the new combined estimator is smaller than that of each individual estimator in and .

The goal is to find weights such that the variance of the new estimator is smaller than the minimum of variance of the estimators by [13] and [17], which have good theoretical properties. To obtain the estimator that yields smallest variance with these properties, we can use the idea of [25], which was applied to the problem of modeling multivariate failure times.

In [25], the joint distribution of estimators is considered, where m = 1,…, k and r = 1…R. In this case, m indicates index of regression parameters and r stands for index of the rth event. For obtaining an optimal estimator, they applied arguments from [24] which derived a linear combination of test statistic to maximize power against every alternative hypothesis. Let be the covariance matrix for the estimators . Then we fix m and define be covariance matrix of . It can be obtained from the entire covariance matrix by selecting the part corresponding to for r = 1,…, R under fixed m. Now we can define , where d = (d 1, d 2,…, d ) satisfies [25]. Then is a vector of weights which leads the best estimator among linear combinations of estimators of where e is a vector consisting of R ones [24], [25].

We now apply the argument in previous paragraph to our model by considering the joint distribution of . Let and where are estimating equations for 0. The strong consistency and asymptotic joint distribution of three estimators, described in following theorems, play a crucial role in our methodology.

To prove asymptotic results, several regularity conditions are required. As stated in [7] and [17], define

Let . Define and

From the Appendix in [17], the additional conditions are as follows: In many parts of proofs, we adapt arguments from [13] and [17].

The parameter space 𝓦 is compact, and the true parameter 0 is an interior point of 𝓦.

0 is the only solution of the estimating equation .

E(||Z||2) < ∞, where ||·|| is Euclidean norm and there exists positive constant K such that partial derivatives of F are bounded by K and there exists positive constant K* such that marginal probability distribution of F is bounded by K* almost surely.

cov[(Z 1 − Z 2){T 1(Z 1, Z 2)}1/2] and cov[(Z 1 − Z 2){T 2(Z 1, Z 2)}1/2] are positive definite.

Theorem 1. By conditions of C1 − C3 in Appendix of [17] and conditions in [27], is (strongly) consistent.

Proof. Let . It suffices to show that and are strongly consistent, respectively. Let = ( , ). Note that we have compact region, say 𝓦 and we assume regularity conditions in [27]. By [27], there exists nonrandom function m 1 such that sup||n −1/2 S () − m 1()|| converges to 0 with probability 1 where 𝓝0 is a neighborhood of 0. Thus is strongly consistent. Similarly, we have another nonrandom function m 2 such that converges to 0 with probability 1 where 𝓝1 is a neighborhood of 0. Hence by [27], is strongly consistent.

For , by argument in Appendix of [17], note that by the U-statistics version of the law of large numbers, for all ∈ 𝓦, converges to 0 in probability where . We can partition our compact space as 𝓦1,…,𝓦 so that . Clearly, then for { ∈ 𝓦, j = 1,…, k}, converges to 0 in probability. Then by Appendix of [17], and for all ϵ > 0, there exists ξ > 0 such that Hence

Thus is strongly consistent and clearly, is strongly consistent.

Theorem 2. Assuming certain technical conditions from [27] and [17], is asymptotically normal with mean zero vector and covariance matrix Σ 0 where where Γ 0 is a nonsingular matrix and Ω 0 is the asymptotic covariance matrix of G ( 0).

Proof. As consistency, we assume the same regularity conditions as in [27]. Let and . Similar to [13], let be the cause-specific hazard function for the and let be the cause-specific hazard function for under dependent censoring. Define Then M 1 and M 2 are martingales [5], [13]. By adapting a proof in the Appendix in [13], Rebdolledo’s martingale central limit theorem [5] gives where and . From Appendix of [17], where 2h 1(v, 0) = 2E[h(v,V 2, 0)]. For j = 1,…, n, M 1(t) is the martingale associated with , while M 2(t) is the martingale associated with and h(V ,V ,) = (Z − Z )ϕ () [13], [17]. For j = 1,…, n, define By the Cramér-Wold theorem, G ( 0) has an asymptotically normal distribution with mean zero and covariance matrix Ω 0, where Note that . As stated in the Appendix of [17], under conditions of N1 − N3 from [9], there exists an open neighborhood of 0, say K 0, such that Using a Taylor series expansion of γ() around 0, With these two results (Eq 5) and (Eq 6), by Appendix of [17], From [27], we have that for any in the small neighborhood of 0, where P 0 is k × k nonsingular matrix. From the Appendix in [13], for , for any in the small neighborhood of 0, where L 10 is defined as is 2k × 2k nonsingular matrix and M 0 and H 0 are k × k constant matrices. Define . Using expansion from [17], for any in the small neighborhood of 0, where and . Combining expansions of (Eq 8), (Eq 9) and (Eq 10), we have for any in the small neighborhood of 0, where Γ 0 is defined as The results from [9] and [27], along with the consistency of , imply that By combining the above results with Slutsky’s theorem, has an asymptotically normal distribution with mean zero and covariance matrix .

Theorem 2 implies the asymptotic normality of with the form of Σ 0 being Let be the estimated covariance matrix of Σ 0. In this covariance matrix, is a k × k covariance matrix for , is a k × k covariance matrix for and is a k × k covariance matrix for . Moreover, and represent covariance terms between and and between and , respectively. Define as the covariance matrix between and . Clearly, , and .

The issue remains of how to obtain the matrix corresponding to in our context. Note that and are correlated with each other. The estimating equation structure implies that and cannot be estimated separately from . Thus our matrix corresponding to should include the effect of . To obtain the matrix, we need to invert whole matrix and extract submatrix corresponding to and . There are two approaches to obtain submatrix.

The first approach is to invert and obtain the submatrix of corresponding to and . Let us denote this matrix as . Clearly, this matrix is 2 × 2 and also positive definite. Then we can calculate , where h = (1, 1). By using the form of the optimal estimator in [25], we obtain new weighted estimator for mth covariate, say , where We can repeat this step for the other regression coefficients. Then we obtain . In this first approach, weights are generated through using k number of 2 × 2 matrices. We can refer this first approach as ‘marginal approach’.

Sometimes it is desirable to consider entire covariates all at once when obtaining weights. The second approach is to obtain the corresponding submatrix of for . We denote this matrix as . This approach is different from first one in that consists of elements of the covariance matrix from and but now has elements of covariance matrix from corresponding entire . This approach reflects the effect of jointly on our new estimator. Let E be a 2k × k matrix such that

E is a multivariate extension of h. Note that E is concatenation of two k × k identity matrices by row. Entries that are 1 in these two k × k identity matrices are source of weights for and . The next step is to construct , which is

Then has the form This matrix is a multivariate extension of from the first approach. This matrix is a contrast matrix in the sense that for the mth regression coefficient of and . Moreover, for p ≠ m = 1,…, k. Using a vector form, from this approach our new estimator, say , We can also refer this approach as the ‘joint approach’.

Now the key step is to obtain . We use the resampling approach of [16], which was also used in [13] and [17]. Let and . From [13] and [17], we have and Define A consistent estimator of Ω 0 is We then solve the estimating equation where Q (i = 1,…, n) represent standard normal random variables. Note that is joint estimating equation for . By solving this equation, we obtain many realizations of s, say where are solutions from (Eq 11). The next theorem, combined with Theorem 2, justifies the resampling approach for calculating .

Theorem 3. Based on the technical conditions in [16], the unconditional distribution of is same asymptotically as the conditional distribution of where are realizations of from resampling.

Proof. Recall that for any in the small neighborhood of 0, we have Note that are solutions of Eq (11). By conditioning on observed data and using expansion (Eq 12) as well as by adapting arguments in [13] and [16], and hence, Note that is asymptotically normal with covariance matrix Σ 0. Then given observed data, distribution of is asymptotically normal with covariance matrix . Hence conditional distribution of on observed data is asymptotically same as unconditional distribution of .

For m = 1,…k and j = 1,…, M, let and be jth realizations of an element and corresponding to mth covariate, respectively. The algorithm for the first approach is as follows.

By resampling, calculate the covariance matrix using realizations and , (m = 1,…, k and j = 1,…, M).

From , obtain the covariance matrix corresponding to and , say .

Calculate where h = (1,1) and obtain the new estimate .

Repeat step 3 for all covariates.

The algorithm for the second approach is as follows.

By resampling, calculate the covariance matrix using realizations and (m = 1,…, k and j = 1,…, M).

Obtain from .

From and E, obtain .

Calculate the new estimate , where be the element of jth row and lth column of .

By Theorem 1 and Theorem 2, our new estimators are consistent and asymptotically normal.

Model checking

For assessing the adequacy of the model, since our weight estimator is based on estimators from [13] and [17], it is reasonable to consider entire processes from [13] and [17]. In this case, we extend model checking technique from [13]. As defined in [13], Let and , where i = 1,…, n. Then Nelson-Aalen estimators for the event of interest and dependent censoring are Note that by (Eq 3) and (Eq 4), martingale residuals are defined as where can be either or . Then as defined in [13],

Then similar to [13] and [17], we can substitute on S (s; ), and on U (t;). are called observed score processes with respect to dependent censoring and the event of interest, respectively [7], [13], [17]. We can construct [13], [17], where where and . These three processes are called bootstrapped processes [7], [13], [17]. We can plot the observed process with bootstrapped processes by randomly selecting 20 or 30 observations. Standard tests for goodness of fit can be performed by calculating Kolmogorov-Smirnov type test statistics. Test statistics are then defined by , and . To calculate the null distribution of the test statistics, first we obtain jth realizations of bootstrap samples and . Then we compute and , respectively for j = 1,…, M, where and are jth realizations of bootstrap samples of and . The p-values can be defined by

[10]. If the p-value is smaller than predetermined level, we reject the null hypothesis, which means that data does not have appropriate fit on our bivariate model. Note that a multiple testing problem arises for testing the models for . We address this by adjusting p-values based on a Bonferroni correction with two tests.

Simulation Studies

We consider two simulation settings. In first simulation setting, the errors follow a bivariate normal distribution with mean (0,1.2) with variance 1 and correlation ρ = 0,0.25. The independent censoring time C is generated from log(U*), where U* has uniform distribution with minimum value 0 and maximum value 20. Covariate is Z ∼ Bernoulli(0.5), where Bernoulli(0.5) is Bernoulli distribution with success probability 0.5. We run 500 simulation runs. Within each simulation run, 500 resampling runs are tried for covariance matrix calculation. Sample sizes are N = 150 and N = 300. If there is only one covariate in the model, the first and the second method of the weighted estimation are equivalent. Let this common weighed estimator be . We calculate bias (Bias), mean squared error (MSE), mean of standard error (SEE), 95% coverage rate (Coverage). The coverage is based on the normal approximation. Moreover, to evaluate robustness of estimators, we also compute median of difference of the estimator from true value (Dmedian), median of squared error of estimates (Mediansq), and median of standard errors (Sdmedian). Results are summarized on Table 1 and Table 2.

Table 1

Simulation result when N = 150 and N = 300, ρ = 0 with covariate Bernoulli(0.5).

Estimators 1	N = 150
	Bias (Dmedian 2 )	MSE 3 (Mediansq 4 )	SEE 5 (Sdmedian 6 )	Coverage 7
θ^^L	0.018 (0.018)	0.04 (0.017)	0.204 (0.2)	0.95
θ^^P	0.021 (0.014)	0.036 (0.017)	0.193 (0.19)	0.96
θ^^WE	0.016 (0.006)	0.036 (0.015)	0.188 (0.185)	0.95
Estimators 1	N = 300
	Bias (Dmedian 2 )	MSE 3 (Mediansq 4 )	SEE 5 (Sdmedian 6 )	Coverage 7
θ^^L	-0.002 (-0.003)	0.017 (0.006)	0.140 (0.140)	0.95
θ^^P	-0.001 (0.002)	0.016 (0.007)	0.133 (0.132)	0.95
θ^^WE	-0.004 (-0.002)	0.016 (0.007)	0.130 (0.129)	0.94

1 : the estimator by [13]; : the estimator by [17]; : the weighted estimator by the proposed approach (Note that the marginal approach and the joint approach are equal in one variable case)

2 median of difference of the estimator from true value

3 mean squared error

4 median of squared error

5 mean of standard error

6 median of standard error

7 95% coverage rate

Table 2

Simulation result when N = 150 and N = 300, ρ = 0.25 with covariate Bernoulli(0.5).

Estimators 1	N = 150
	Bias (Dmedian 2 )	MSE 3 (Mediansq 4 )	SEE 5 (Sdmedian 6 )	Coverage 7
θ^^L	0.005 (0.01)	0.036 (0.017)	0.198 (0.197)	0.95
θ^^P	0.006 (0.007)	0.032 (0.015)	0.189 (0.188)	0.95
θ^^WE	-0.001 (-0.006)	0.033 (0.016)	0.184 (0.183)	0.94
Estimators	N = 300
	Bias (Dmedian 2 )	MSE 3 (Mediansq 4 )	SEE 5 (Sdmedian 6 )	Coverage 7
θ^^L	-0.003 (0.005)	0.018 (0.008)	0.138 (0.137)	0.95
θ^^P	0.001 (0.007)	0.017 (0.007)	0.131 (0.131)	0.95
θ^^WE	-0.003 (0.002)	0.017 (0.007)	0.129 (0.128)	0.95

1 : the estimator by [13]; : the estimator by [17]; : the weighted estimator by the proposed approach (Note that the marginal approach and the joint approach are equal in one variable case)

2 median of difference of the estimator from true value

3 mean squared error

4 median of squared error

5 mean of standard error

6 median of standard error

7 95% coverage rate

In second simulation setting, we generate Gamma random variable ν with mean μ = 1 and variance σ 2 = 0 or 1, then create W = exp (ϵ ), which is an exponential random variable with rate 4ν −1 and exp (ϵ ) with an exponential random variable with rate ν −1. Then we generate time to the event of interest by and time to the dependent censoring by (By notation in our paper, X, D and C are already log-transformed times. Thus in this context, exp (X), exp (D) and exp (C) are times in the original scale). The independent censoring time exp (C) has uniform distribution with minimum value 0 and maximum value 20. True parameter values are 0 = (0.5,1) and 0 = (1,0.5) and covariates are Z 1 ∼ U(0,1), where U(0,1) is uniform distribution with minimum value 0 and maximum value 1 and Z 2 ∼ Bernoulli(0.5). We run 500 simulation runs. Within each simulation run, 500 resampling runs are tried for covariance matrix calculation. Let be weighted estimators from calculating weights marginally (the first proposed method) and let be weighted estimators from calculating weights jointly (the second proposed method). We compute the same quantities as we did in the first set of the simulation study. Results are summarized on Table 3 and Table 4.

Table 3

Simulation result when N = 150 and N = 300, σ2 = 0 with two covariates (Z1: U(0, 1), Z2: Bernoulli(0.5)).

Estimators 1	N = 150
	Bias (Dmedian 2 )		MSE 3 (Mediansq 4 )		SEE 5 (Sdmedian 6 )		Coverage 7
	Z 1	Z 2	Z 1	Z 2	Z 1	Z 2	Z 1	Z 2
θ^^L	0.0001 (0.002)	0.002 (-0.005)	0.12 (0.052)	0.042 (0.018)	0.358 (0.352)	0.226 (0.222)	0.96	0.96
θ^^P	-0.003 (0.003)	-0.003 (-0.002)	0.158 (0.074)	0.051 (0.023)	0.427 (0.418)	0.243 (0.241)	0.96	0.96
θ^^MWE	0.003 (-0.007)	0.003 (0.003)	0.123 (0.053)	0.043 (0.019)	0.351 (0.349)	0.219 (0.218)	0.95	0.95
θ^^JWE	0.003 (-0.007)	0.004 (0.001)	0.123 (0.055)	0.043 (0.018)	0.351 (0.349)	0.219 (0.217)	0.94	0.95
Estimators 1	N = 300
	Bias (Dmedian 2 )		MSE 3 (Mediansq 4 )		SEE 5 (Sdmedian 6 )		Coverage 7
	Z 1	Z 2	Z 1	Z 2	Z 1	Z 2	Z 1	Z 2
θ^^L	-0.012 (-0.013)	0.004 (0.001)	0.065 (0.028)	0.02 (0.01)	0.257 (0.255)	0.148 (0.148)	0.95	0.96
θ^^P	-0.014 (-0.017)	-0.001 (-0.015)	0.081 (0.035)	0.024 (0.013)	0.283 (0.281)	0.162 (0.162)	0.95	0.96
θ^^MWE	-0.01 (-0.012)	0.003 (-0.001)	0.064 (0.031)	0.02 (0.01)	0.252 (0.251)	0.146 (0.146)	0.95	0.96
θ^^JWE	-0.01 (-0.014)	0.003 (0.002)	0.064 (0.032)	0.02 (0.01)	0.251 (0.25)	0.146 (0.146)	0.95	0.96

1 : the estimator by [13]; : the estimator by [17]; : the weighted estimator by the marginal approach; : the weighted estimator by the joint approach

2 median of difference of the estimator from true value

3 mean squared error

4 median of squared error

5 mean of standard error

6 median of standard error

7 95% coverage rate

Table 4

Simulation result when N = 150 and N = 300, σ2 = 1 with two covariates (Z1: U(0,1), Z2: Bernoulli(0.5)).

Estimators 1	N = 150
	Bias (Dmedian 2 )		MSE 3 (Mediansq 4 )		SEE 5 (Sdmedian 6 )		Coverage 7
	Z 1	Z 2	Z 1	Z 2	Z 1	Z 2	Z 1	Z 2
θ^^L	-0.010 (-0.003)	-0.033 (-0.039)	0.273 (0.127)	0.086 (0.038)	0.441 (0.434)	0.315 (0.312)	0.90	0.96
θ^^P	0.002 (-0.002)	-0.032 (-0.046)	0.295 (0.142)	0.095 (0.038)	0.559 (0.552)	0.325 (0.324)	0.95	0.96
θ^^MWE	-0.009 (-0.008)	-0.030 (-0.031)	0.263 (0.128)	0.085 (0.041)	0.437 (0.432)	0.303 (0.301)	0.90	0.96
θ^^JWE	-0.009 (-0.008)	-0.030 (-0.03)	0.262 (0.128)	0.086 (0.04)	0.437 (0.432)	0.303 (0.301)	0.90	0.96
Estimators 1	N = 300
	Bias (Dmedian 2 )		MSE 3 (Mediansq 4 )		SEE 5 (Sdmedian 6 )		Coverage 7
	Z 1	Z 2	Z 1	Z 2	Z 1	Z 2	Z 1	Z 2
θ^^L	-0.024 (-0.038)	0.003 (-0.006)	0.133 (0.057)	0.04 (0.019)	0.345 (0.344)	0.211 (0.21)	0.94	0.96
θ^^P	-0.016 (-0.016)	0.012 (0.003)	0.148 (0.059)	0.045 (0.02)	0.384 (0.382)	0.222 (0.222)	0.96	0.97
θ^^MWE	-0.024 (-0.035)	0.007 (-0.003)	0.134 (0.058)	0.04 (0.019)	0.341 (0.341)	0.207 (0.207)	0.94	0.96
θ^^JWE	-0.025 (-0.035)	0.007 (-0.002)	0.135 (0.058)	0.039 (0.018)	0.341 (0.341)	0.206 (0.207)	0.94	0.96

1 : the estimator by [13]; : the estimator by [17]; : the weighted estimator by the marginal approach; : the weighted estimator by the joint approach

2 median of difference of the estimator from true value

3 mean squared error

4 median of squared error

5 mean of standard error

6 median of standard error

7 95% coverage rate

In these simulation results, we can see that our weighted estimators have good results. In both cases, bias and mean squared error of our new estimator has similar performance compared to the estimators by [13] and [17]. Mean of standard errors and median of standard errors are smaller than the estimators by [13] and [17]. Moreover, computation results for the median of difference of the estimators from true value and the median of squared error imply that our proposed estimator is comparable with the estimators from the original methods.

In the first simulation setting, the difference of standard error between our proposed estimator and is bigger than the one between and the proposed estimator. In the second simulation setting, the phenomenon is the opposite. Furthermore, in the first simulation setting, has lower standard error on average than one of while have better efficiency (with respect to standard error) than ones by in the second simulation setting. This simulation result verifies our claim, which means that neither estimator is better than another. Our proposed estimator takes advantage of smaller standard error with achieving small bias and correct coverage except N = 150 with σ 2 = 1 in the second simulation setting. In this scenario, empirical coverage of proposed estimators is lower than nominal 95% coverage. This is due to low coverage of . Since we combine and , if one of them has low coverage, it is highly likely that the coverage of weighted estimator may also be below the nominal coverage.

Real data analysis

We applied our method to data from the AIDS Clinical Trial Group (ACTG) Study 364 [1], which was used in [17]. This multicenter randomized study investigated patients whose plasma RNA level is at least 500 copies per ml. Subjects were assigned to three treatments, nelfinavir (NFV), efavirenz (EFV), and combination of nelfinavir and efavirenz (NFV + EFV). Details about this study can be found in [1].

The two failure times are time to HIV RNA level greater than 2000 copies per ml and time to withdrawal of study. Let X be the first time when HIV RNA level is greater than 2000 copies per ml and D be time to withdrawal of study. We considered four covariates and 194 observations. Z 1 takes value 1 if a patient receives EFV and 0 otherwise. Z 2 takes value 1 if a patient receives NFV + EFV and 0 otherwise. Z 3 is New3TC, which takes value 1 if lamivudine is given as a new nucleoside analogue therapy to a patient and 0 otherwise. Z 4 is logarithm of RNA level at the start of the study.

Table 5 and Table 6 show the point estimates and standard errors of , , , and . Our method works well for the models with and without New3TC on all covariates. Some variables are seen to be statistically significant based on the weighted estimator while they are not by [13] or [17]. For example, let’s consider effect of EFV to the time to first virologic failure. By Table 6, the estimated effect by using approach of [13] is 0.475 and its standard error is 0.250. From approach of [17], an estimate is 0.464 and its standard error is 0.281. Based on the fact that estimators are asymptotic normal, from Wald test using [13] and [17], EFV is not a statistically significant variable on 5% significant level. On the other hand, a weighted estimate using first approach is 0.471 and its standard error is 0.222. In this case, EFV is a statistically significant variable on 5% significant level.

Table 5

Point estimates with standard errors of covariates in AIDS study for model without New3TC (Standard errors are shown in parenthesis).

Covariates	η^	θ^^L	θ^^P	θ^^MWE	θ^^JWE
EFV 1	0.753 (0.339)	0.115 (0.219)	0.375 (0.269)	0.168 (0.206)	0.2 (0.205)
NFV 2 + EFV	0.674 (0.255)	1.128 (0.239)	1.091 (0.309)	1.120 (0.222)	1.114 (0.222)
log(RNA) 3	-0.544 (0.154)	-0.464 (0.215)	-0.531 (0.169)	-0.507 (0.163)	-0.511 (0.162)

1 efavirenz

2 nelfinavir

3 logarithm of RNA at the start of the study

Table 6

Point estimates with standard errors of covariates in AIDS study for model with New3TC (Standard errors are shown in parenthesis).

Covariates	η^	θ^^L	θ^^P	θ^^MWE	θ^^JWE
EFV 1	0.770 (0.278)	0.475 (0.250)	0.464 (0.281)	0.471 (0.222)	0.471 (0.222)
NFV 2 + EFV	0.650 (0.260)	1.353 (0.277)	1.246 (0.338)	1.333 (0.263)	1.317 (0.261)
New3TC 3	0.927 (0.355)	1.449 (0.296)	1.374 (0.328)	1.431 (0.267)	1.420 (0.261)
log(RNA) 4	-0.631 (0.183)	-0.654 (0.289)	-0.661 (0.218)	-0.659 (0.216)	-0.660 (0.215)

1 efavirenz

2 nelfinavir

3 lamivudine as new nucleoside analogue therapy

4 logarithm of RNA at the start of the study

Observed score process with bootstrapped processes for withdrawal of study with respect to Z 1 is shown in Fig 1. Fig 2 and Fig 3 show observed score processes and bootstrapped processes of the first virologic failure using with respect to Z 1. These three plots are based on the model without New3TC. They are fluctuating around zero, so it seems that there is no graphical evidence for lack of fit. The p-value for the lack of fit tests of withdrawal is 0.952 and the first virologic failure using and are 0.918 and 0.959 respectively. With graphical checking, p-value indicates that there is no evidence for violation of the model assumption.

Fig 1

Plot of observed score process and bootstrapped processes of time to withdrawal of study with respect to Z 1.

The thickline is observed process and the dashed lines are bootstrapped processes.

Fig 2

Plot of observed score process and bootstrapped processes of time to first virologic failure using with respect to Z 1.

The thickline is observed process and the dashed lines are bootstrapped processes.

Fig 3

Plot of observed score process and bootstrapped processes of time to first virologic failure using with respect to Z 1.

The thickline is observed process and the dashed lines are bootstrapped processes.

For purposes of interpretation, since D represents a standard survival time, the interpretation of is in terms of covariate effect for survival time. However, since the observed time for X depends on D, interpretation of is difficult. One way to interpret is to assume that D does not exist and interpret the effect of on X only. This approach is possible if there exists a reasonable extrapolation mechanism for X [18]. However, considering the estimation structure for , it is difficult to separate effect of to X from effect of to D.

Discussion

In this paper, we have proposed optimal estimators using combinations of the two estimators from [13] and [17]. Our methodology can be extended to a case of recurrent event with dependent censoring, which is extensively studied [6], [7], [10]. We are currently working on this extension.

Optimality of the estimator has been discussed in other contexts. Recently, there is a publication that proposed optimal additive functions based on score functions [14]. The main point of their method is to combine unbiased estimating functions. In our case, this would be combining estimating equations and new solution can be obtained by this estimating equation. Comparing performance of this solution and our proposed estimator is of interest. This will be left open to future research.

Another way of achieving optimality is to use generalized method of moment estimator [8]. This estimator is a linear combination of estimating functions [19]. In this case, the estimating functions have a greater dimension than the dimension of the parameter vector. The optimality is achieved by the linear combination. It is shown that the estimator from this linear combination of estimating functions is consistent and asymptotically normal [8]. In the literature of statistics, this idea is applied to generalized estimating equations [19]. The estimating functions proposed by [19] are called quadratic inference function. Recently, the quadratic inference function is applied to Cox model [26].

[8] and [19] derived new estimating functions, while we combined two estimators directly. This idea of the generalized method of moments is very appealing, but the estimating functions of [13] and [17] are nonsmooth. Finding derivative for the linear combination of the estimating functions, which is a key in generalized method moments, is challenging for our work because we cannot find the derivatives in the estimating functions proposed by [13] and [17]. Applying the idea of [8] to AFT model will be interesting future research.

Our estimating equations to obtain estimators involve nonsmooth functions of and . Many literatures used a linear programming approach for estimating [3], [11]. However, this linear programming method is very slow for computing estimators of . Thus this approach is very inefficient when implementing to solve (Eq 11) for estimation of Σ. Recently, an approach called a derivative free-spectral algorithm for nonlinear equations (DF-SANE) was proposed [12], and there is a publication that showed that this algorithm is better than the linear programming method using an example of estimating parameters of AFT models under independent censoring. [21]. However, under dependent censoring, the artificial censoring term leads to numerical instability in estimating parameters and calculating resampled estimators. Moreover, this algorithm does not converge well under default tolerance settings using DF-SANE [21]. Thus using this algorithm requires changing the tolerance level. Developing efficient numerical algorithms for estimating parameters is an important topic for future research.