# An Efficient Estimation of the Mean Residual Life Function with Length-Biased Right-Censored Data

Authors: Hongping Wu; Yihui Luan
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The mean residual life (MRL) function for a lifetime random variableT0 is one of the basic parameters of interest in survival analysis. In this paper, we propose a new estimator of the MRL function with length-biased right-censored data and evaluate its performance through a small Monte Carlo simulation study. The results of the simulations show that the proposed estimator outperforms the existing one referred to in Data and Model Setup Section in terms of Monte Carlo bias and mean square error, especially when the censoring rate is heavy. We also show that the proposed estimator converges in distribution under some conditions.


---

## Body

## 1. Introduction

The MRL function at timet is defined to be the expected remaining lifetime of a system given survival up to time t. It can be used to model various lifetime data in many areas of science. For example, a life insurance company may be interested in the life expectancy of a person, or an engineering firm may wish to estimate the expected remaining lifetime of a system, given survival up to time t. In literature, there has been a lot of work on the inference of the MLR function. Based on the i.i.d copies of the life time and the closed relationship between the MRL function and survival function, nonparametric estimators and NA based confidence intervals have been constructed (see [1, 2]). For right-censored data, consistent estimator of the MRL function with its asymptotic normal distribution has been given by Yang [3] and Kumazawa [4]. EL approach has been developed to construct the confidence interval of the MRL function by Zhao and Qin [5], where, for censored data, the corresponding log-EL ratio was defined and the limiting distribution was proved to be a scaled χ2 distribution. A class of transformed MRL function for fitting survival data under right censoring was proposed and efficiency and robustness of estimators have been studied by Sun and Zhang [6]. A more recent work on estimation of the MRL with left-truncated and right-censored data has been constructed by Zhao et al. [7], and they showed that the proposed estimator converges weakly to a Gaussian process.This paper considers survival data arising from length-biased sampling, which is a special case of left truncation. Following the terminology in the literature (see [8, 9]), length-biased data are defined for left-truncated and right-censored data under the stationarity assumption, which means that the probability of the survival time being sampled is proportional to its length or the survival times are left-truncated by uniformly distributed random truncation times. Statistical inferences of length-biased data have been considered by many authors. Vardi [9] considered the nonparametric maximum likelihood estimator of the survival distribution. Gupta and Keating [10] obtained the relations for reliability measures of the length-biased distribution and some characterization results. Concerning the effect of the length bias on the sampling distribution of the covariates, Bergeron et al. [11] presented a joint likelihood approach and studied the large-sample behavior of the resulting maximum likelihood estimators. To estimate the survival distribution, Luo and Tsai [12] proposed a pseudopartial likelihood approach, and Huang and Qin [8] proposed a new nonparametric estimator based on the product-limit estimator. Their estimators were proved to be uniformly consistent and converge weakly to Gaussian processes.In this paper we propose a new estimator of the MRL function with length-biased and right-censored data. We compare the performance of our estimator with the existing one referred to in Section2.1 through a small Monte Carlo simulation study and the results of the simulations show that the proposed estimator outperforms the existing one in terms of Monte Carlo bias and mean square error, especially when the censoring rate is heavy. The research is organized as follows. In Section 2, we introduce the model and the proposed estimator of the MRL function. The large sample property of our proposed estimator is stated in Section 3. The performance of the estimators is illustrated in Section 4.


## 2. Notations and Estimators

### 2.1. Data and Model Setup Section

LetA0 and T0 be two independent nonnegative random variables denoting the left-truncated time and the lifetime of a system, respectively. Assume that the truncation probability α=Pr(A0≤T0)>0. Let A and T denote the observed truncation time and lifetime; that is, the pair (A,T) is a realization from the conditional distribution of (A0,T0) given that A0≤T0. Under the stationarity assumption, A is uniformly distributed on [0,T], so T has a length-biased distribution. Coupled with the censoring time C, which is applied to the residual lifetime V=T-A instead of T and independent of (A,T), the observed data is composed of i.i.d. copies of (A,Y,δ), where Y=min⁡(T,A+C), δ=I(V≤C). Now we suppose {(ai,yi,δi),i=1,2,…,n} to be a sequence of i.i.d. random samples of (A,Y,δ) which one observes.Denote byf(t) and S(t) the marginal density function and survival function of T0 and by g(t) and G(t) the corresponding density function and survival function of the length-biased random variable T. Let τ=inf⁡{t:S(t)=0}, τ∈R+∪{+∞}. Under length-biased sampling, A and T are stochastically dependent and (A,T) has the joint density function (see [13])
(1)fA,T(a,t)=f(t)μ,0≤a<t<τ,
and the density function f of T0 is related to the length-biased density function g of T, as follows:
(2)g(t)=tf(t)μ,0≤t<τ,
where μ=[∫0τ‍u-1g(u)du)]-1 is finite.According to the definition of the MRL function, the mean residual life ofT0 at time 0≤t<τ is expressed as
(3)m(t)=E(T0-t∣T0>t)=∫tτ‍S(u)duS(t),
and m(t)=0 when t≥τ.In Huang and Qin [8], the nonparametric maximum likelihood estimator of S(t) has been given:
(4)S^(t)=∏u∈[0,t]‍{1-dΛ^(u)},
where Λ^ is the estimated cumulative hazard function
(5)Λ^(t)=∫0t‍dN¯(u)R¯(u),
with N¯(u)=n-1∑i=1n‍δiI(yi≤u), R¯(u)=n-1∑i=1n‍I(yi≥u)-S^A(u), S^A(u)=∏v∈[0,u]‍{1-(dQ~(v)/K~(v))}, Q~(v)=n-1∑i=1n‍I{I(ai≤v)+δiI(yi-ai≤v)}, and K~(v)=n-1∑i=1n‍{I(ai≥v)+I(yi-ai≥v)}.Applying (3) and (4), when 0≤t≤y(n), the estimator of m(t) can be naturally derived:
(6)m^(t)=∫ty(n)‍S^(u)duS^(t),
where y(n)=max⁡(y1,y2,…,yn). If t>y(n), m^(t)=0.LetSC(·) be the survival function of C and τC=inf⁡{t:SC(t)=0}. Assume the following conditions hold:(i)
L is a nonnegative constant satisfying L<sup⁡{t:S(t)SC(t)>0};(ii)
S and SC are absolutely continuous on [0,L];(iii)
∫tτ‍g(u)du<∞, ∫tτ‍u-1g(u)du<∞, t∈[0,τ);(iv)
eitherτC<τ with SC(τC-)>0 or τC≥τ.As a straightforward result of the functional delta method, whenn→∞, n(m^(t)-m(t)) converges in distribution to a zero mean normal variable with the variance function σ2(t)=S-4(t)[S2(t)E(∫tτ‍Z(u)du)2-2S(t)∫tτ‍Z(u)duEZ(t)∫tτ‍Z(u)du+(∫tτ‍S(u)du)2EZ2(t)] for every t∈[0,L] under Assumptions (i–iv).In the next subsection, we shall state our proposed estimator ofm(t), and the simulation study in Section 4 shows that it performs better in several terms than the one in (6).


### 2.2.  The Proposed Estimator

According to the joint density function of(A,T), we can derive that
(7)E[I(Y≥t)δSC(Y-A)]=∫tτ‍f(u)μ∫0u‍(∫u-a∞‍fC(v)dv)1SC(u-a)dadu=∫tτ‍g(u)du,E[I(Y≥t)δYSC(Y-A)]=∫tτ‍f(u)μ∫0u‍(∫u-a∞‍fC(v)dv)1u1SC(u-a)dadu=∫tτ‍u-1g(u)du,
where fC(·) is the density function of the censoring time C.Applying (7) we shall use the method of moments to derive an estimator for the MRL function m(t). To proceed we note that m(t) can be written as
(8)m(t)=E[T0-t∣T0>t]=∫tτ‍uf(u)du∫tτ‍f(u)du-t=∫tτ‍g(u)du∫tτ‍u-1g(u)du-t.An alternative estimator of the MLR functionm(t), for 0≤t≤y(n), can be derived as follows:
(9)m~(t)=(1/n)∑i=1n‍I(yi>t)δi/S^C(yi-ai)(1/n)∑i=1n‍yi-1I(yi>t)δi/S^C(yi-ai)-t,
where S^C(t) is the product-limit estimation of SC(t), and it is expressed as
(10)S^C(t)=∏i=1n‍(1-(1-δi)I(yi-ai≤t)R(yi-ai)),
where R(t)=∑j=1n‍I(yj-aj≥t). If t>y(n), m~(t)=0.In the next section we shall establish thatm~(t) converges in distribution under some conditions by using the Hadamard derivative method.


## 2.1. Data and Model Setup Section

LetA0 and T0 be two independent nonnegative random variables denoting the left-truncated time and the lifetime of a system, respectively. Assume that the truncation probability α=Pr(A0≤T0)>0. Let A and T denote the observed truncation time and lifetime; that is, the pair (A,T) is a realization from the conditional distribution of (A0,T0) given that A0≤T0. Under the stationarity assumption, A is uniformly distributed on [0,T], so T has a length-biased distribution. Coupled with the censoring time C, which is applied to the residual lifetime V=T-A instead of T and independent of (A,T), the observed data is composed of i.i.d. copies of (A,Y,δ), where Y=min⁡(T,A+C), δ=I(V≤C). Now we suppose {(ai,yi,δi),i=1,2,…,n} to be a sequence of i.i.d. random samples of (A,Y,δ) which one observes.Denote byf(t) and S(t) the marginal density function and survival function of T0 and by g(t) and G(t) the corresponding density function and survival function of the length-biased random variable T. Let τ=inf⁡{t:S(t)=0}, τ∈R+∪{+∞}. Under length-biased sampling, A and T are stochastically dependent and (A,T) has the joint density function (see [13])
(1)fA,T(a,t)=f(t)μ,0≤a<t<τ,
and the density function f of T0 is related to the length-biased density function g of T, as follows:
(2)g(t)=tf(t)μ,0≤t<τ,
where μ=[∫0τ‍u-1g(u)du)]-1 is finite.According to the definition of the MRL function, the mean residual life ofT0 at time 0≤t<τ is expressed as
(3)m(t)=E(T0-t∣T0>t)=∫tτ‍S(u)duS(t),
and m(t)=0 when t≥τ.In Huang and Qin [8], the nonparametric maximum likelihood estimator of S(t) has been given:
(4)S^(t)=∏u∈[0,t]‍{1-dΛ^(u)},
where Λ^ is the estimated cumulative hazard function
(5)Λ^(t)=∫0t‍dN¯(u)R¯(u),
with N¯(u)=n-1∑i=1n‍δiI(yi≤u), R¯(u)=n-1∑i=1n‍I(yi≥u)-S^A(u), S^A(u)=∏v∈[0,u]‍{1-(dQ~(v)/K~(v))}, Q~(v)=n-1∑i=1n‍I{I(ai≤v)+δiI(yi-ai≤v)}, and K~(v)=n-1∑i=1n‍{I(ai≥v)+I(yi-ai≥v)}.Applying (3) and (4), when 0≤t≤y(n), the estimator of m(t) can be naturally derived:
(6)m^(t)=∫ty(n)‍S^(u)duS^(t),
where y(n)=max⁡(y1,y2,…,yn). If t>y(n), m^(t)=0.LetSC(·) be the survival function of C and τC=inf⁡{t:SC(t)=0}. Assume the following conditions hold:(i)
L is a nonnegative constant satisfying L<sup⁡{t:S(t)SC(t)>0};(ii)
S and SC are absolutely continuous on [0,L];(iii)
∫tτ‍g(u)du<∞, ∫tτ‍u-1g(u)du<∞, t∈[0,τ);(iv)
eitherτC<τ with SC(τC-)>0 or τC≥τ.As a straightforward result of the functional delta method, whenn→∞, n(m^(t)-m(t)) converges in distribution to a zero mean normal variable with the variance function σ2(t)=S-4(t)[S2(t)E(∫tτ‍Z(u)du)2-2S(t)∫tτ‍Z(u)duEZ(t)∫tτ‍Z(u)du+(∫tτ‍S(u)du)2EZ2(t)] for every t∈[0,L] under Assumptions (i–iv).In the next subsection, we shall state our proposed estimator ofm(t), and the simulation study in Section 4 shows that it performs better in several terms than the one in (6).


## 2.2.  The Proposed Estimator

According to the joint density function of(A,T), we can derive that
(7)E[I(Y≥t)δSC(Y-A)]=∫tτ‍f(u)μ∫0u‍(∫u-a∞‍fC(v)dv)1SC(u-a)dadu=∫tτ‍g(u)du,E[I(Y≥t)δYSC(Y-A)]=∫tτ‍f(u)μ∫0u‍(∫u-a∞‍fC(v)dv)1u1SC(u-a)dadu=∫tτ‍u-1g(u)du,
where fC(·) is the density function of the censoring time C.Applying (7) we shall use the method of moments to derive an estimator for the MRL function m(t). To proceed we note that m(t) can be written as
(8)m(t)=E[T0-t∣T0>t]=∫tτ‍uf(u)du∫tτ‍f(u)du-t=∫tτ‍g(u)du∫tτ‍u-1g(u)du-t.An alternative estimator of the MLR functionm(t), for 0≤t≤y(n), can be derived as follows:
(9)m~(t)=(1/n)∑i=1n‍I(yi>t)δi/S^C(yi-ai)(1/n)∑i=1n‍yi-1I(yi>t)δi/S^C(yi-ai)-t,
where S^C(t) is the product-limit estimation of SC(t), and it is expressed as
(10)S^C(t)=∏i=1n‍(1-(1-δi)I(yi-ai≤t)R(yi-ai)),
where R(t)=∑j=1n‍I(yj-aj≥t). If t>y(n), m~(t)=0.In the next section we shall establish thatm~(t) converges in distribution under some conditions by using the Hadamard derivative method.


## 3. The Asymptotic Properties

Firstly, let us introduce some notations and lemmas.Denote(11)M(t)=∫tτ‍g(u)du,N(t)=∫tτ‍u-1g(u)du,M^(t)=(1/n)∑i=1n‍I(yi>t)δiS^C(yi-ai),N^(t)=(1/n)∑i=1n‍yi-1I(yi>t)δiS^C(yi-ai),γ(M(t),N(t))=m(t)=M(t)N(t)-t,γ(M^(t),N^(t))=m~(t)=M^(t)N^(t)-t,
and assume that
(12)∬1SC(t-a)F*(da,dt)<∞,∬1t2SC(t-a)F*(da,dt)<∞,
where F*(·,·) is the joint distribution function of (A,T).Lemma1 gives the Hadamard derivative of γ at (M(t),N(t)).Lemma 1.
Fort∈[0,L], the function γ is Hadamard-differentiable at (M(t),N(t)), where its derivative is equal to
(13)γ(M(t),N(t))′(h1(·),h2(·))=h1(t)N(t)-M(t)h2(t)N2(t).Lemma 2.
Assume the conditions (i–iv) and (12) hold, for every t∈[0,L], as n→∞. Then(i)
Z1n(t)=^n(M^(t)-M(t)) converges in distribution to a zero mean normal variable Z1 with the variance function E(I(Y≥t)δ/SC2(Y-A))-M2(t);(ii)
Z2n(t)=^n(N^(t)-N(t)) converges in distribution to a zero mean normal variable Z2 with the variance function E(I(Y>t)δ/Y2SC2(Y-A))-N2(t).Proof.
Firstly, let us give the proof of the result in Lemma2(i). Under the regular conditions (i–iv) and (12), since E(I(Y>t)δ/SC(Y-A))=M(t), by the central limit theorem, we can derive
(14)n(1n∑i=1n‍I(yi>t)δiSC(yi-ai)-M(t))⟶LN(0,EI(Y>t)δSC2(Y-A)-M2(t)),
where →L denotes weak convergence.
Since(15)|1n∑i=1n‍I(yi>t)δiS^C(yi-ai)-M(t)|=|1n∑i=1n‍I(yi>t)δiS^C(yi-ai)-1n∑i=1n‍I(yi>t)δiSC(yi-ai)mmm+1n∑i=1n‍I(yi>t)δiSC(yi-ai)-M(t)|≤|1n∑i=1n‍I(yi>t)δiS^C(yi-ai)-1n∑i=1n‍I(yi>t)δiSC(yi-ai)|+|1n∑i=1n‍I(yi>t)δiSC(yi-ai)-M(t)|≤1n∑i=1n‍I(yi>t)δi|1S^C(yi-ai)-1SC(yi-ai)|+|1n∑i=1n‍I(yi>t)δiSC(yi-ai)-M(t)|,
and recall that S^C(·) is the product-limit estimator of SC(·), so on the basis of (14) and the consistency of S^C(·) (see [14]), we can obtain
(16)|1n∑i=1n‍I(yi>t)δiS^C(yi-ai)-M(t)|⟶0,n⟶∞,Var⁡(Z1n)=EI(yi>t)δiS^C2(yi-ai)-M2(t)⟶EI(Y>t)δSC2(Y-A)-M2(t),n⟶∞.
SinceI(yi>t)δi/S^C(yi-ai), i=1,2,…,n, are i.i.d, it is easy to derive the result of Lemma 2 (i) by the central limit theorem.
Using the same method, let us give the proof of Lemma2 (ii). Under the regular conditions (i–iv) and (12), because E(I(Y>t)δ/YSC(Y-A))=N(t), by the central limit theorem, we can derive
(17)n(1n∑i=1n‍I(yi>t)δiyiSC(yi-ai)-N(t))⟶LN(0,EI(Y>t)δY2SC2(Y-A)-N2(t)).
Since(18)|1n∑i=1n‍I(yi>t)δiyiS^C(yi-ai)-N(t)|=|1n∑i=1n‍I(yi>t)δiyiS^C(yi-ai)-1n∑i=1n‍I(yi>t)δiyiSC(yi-ai)mm+1n∑i=1n‍I(yi>t)δiyiSC(yi-ai)-N(t)|≤|1n∑i=1n‍I(yi>t)δiyiS^C(yi-ai)-1n∑i=1n‍I(yi>t)δiyiSC(yi-ai)|+|1n∑i=1n‍I(yi>t)δiyiSC(yi-ai)-N(t)|≤1n∑i=1n‍I(yi>t)δiyi-1|1S^C(yi-ai)-1SC(yi-ai)|+|1n∑i=1n‍I(yi>t)δiyiSC(yi-ai)-N(t)|,
so on the basis of (17) and the consistency of S^C(·), we can derive
(19)|1n∑i=1n‍I(yi>t)δiyiS^C(yi-ai)-N(t)|⟶0,n⟶∞,Var⁡(Z2n)=EI(yi>t)δiyi2S^C2(yi-ai)-N2(t)⟶EI(Y>t)δY2SC2(Y-A)-N2(t),n⟶∞.
SinceI(yi>t)δi/yiS^C(yi-ai), i=1,2,…,n, are i.i.d, it is easy to derive the result of Lemma 2 (ii) by the central limit theorem. The proof of Lemma 2 is completed.On the basis of Lemmas1 and 2, now we can discuss the asymptotic property of m~(t).Theorem 3.
Assume the conditions (i–iv) and (12) hold, for every t∈[0,L], as n→∞, then n(m~(t)-m(t)) converges in distribution to a zero mean normal variable with the variance function (EZ12(t)/N2(t))-(M(t)EZ1(t)Z2(t)/N3(t))+(M2(t)EZ22(t)/N4(t)).Proof.
Under the conditions (i–iv) and (12), Lemma 2 allows a straightforward application of Lemma 1. So we can derive
(20)n(m~(t)-m(t))⟶dNm~=Z1(t)N(t)-M(t)Z2(t)N2(t),
where →d denotes convergence in distribution. Nm~ is Gaussian as a linear form on two Gaussian variables and it is easy to check that ENm~=0, since Z1 and Z2 are centered. It is also easy to calculate the variance of Nm~ and it can be expressed as
(21)Var⁡(Nm~)=EZ12(t)N2(t)-M(t)EZ1(t)Z2(t)N3(t)+M2(t)EZ22(t)N4(t).
The proof is completed.


## 4. Simulations

To compare the performance of our proposed estimator with the one in (6), we conducted a series of Monte Carlo simulations with n=100 and 200. In these simulations, the underlying survival time T0 has a uniform (0,1) distribution. Then the length-biased density function of T is g(t)=2t, t∈(0,1), and the left-truncated time A is uniformly distributed on [0,T]. We simulated the censoring time C from either a uniform (0,1), a uniform (0,5/3) or a uniform (0,5) distribution, corresponding to a censoring rate of 50%, 30%, and 10%.Table1 summarizes the Monte Carlo bias, standard deviation, and mean square error for each estimator at t=0.2,0.4,0.6 and 0.8 based on 1000 replications. The survival probabilities at the selected time-points are 0.8, 0.6, 0.4, and 0.2. The corresponding true values of m(t) at the selected time-points are 0.4, 0.3, 0.2, and 0.1, respectively.Table 1
Simulation results for nonparametric estimators based on1000 replications.
t
m
(
t
)
m
^
(
t
)
m
~
(
t
)
Bias
SE
MSE
Bias
SE
MSE
n
=
100
,
C
~
Un
(
0,5
)
0.2
0.4
51
28
34
2
29
8
0.4
0.3
50
21
30
0
20
4
0.6
0.2
50
16
28
0
15
2
0.8
0.1
50
14
27
0
11
1
n
=
200
,
C
~
Un
(
0,5
)
0.2
0.4
51
20
30
1
20
4
0.4
0.3
51
14
28
1
14
2
0.6
0.2
50
11
26
0
11
1
0.8
0.1
50
7
25
0
7
1
n
=
100
,
C
~
Un
(
0,5
/
3
)
0.2
0.4
49
32
35
1
30
9
0.4
0.3
49
28
31
1
22
5
0.6
0.2
48
29
31
0
17
3
0.8
0.1
45
56
51
0
12
1
n
=
200
,
C
~
Un
(
0,5
/
3
)
0.2
0.4
50
21
30
1
21
5
0.4
0.3
50
16
28
0
15
2
0.6
0.2
50
15
28
1
12
1
0.8
0.1
51
15
28
1
8
1
n
=
100
,
C
~
Un
(
0,1
)
0.2
0.4
49
62
62
2
32
10
0.4
0.3
53
233
572
−2
26
7
0.6
0.2
37
160
268
5
21
5
0.8
0.1
107
1733
30130
−1
15
2
n
=
200
,
C
~
Un
(
0,1
)
0.2
0.4
49
24
30
2
22
5
0.4
0.3
50
43
43
2
19
4
0.6
0.2
48
26
30
−2
14
2
0.8
0.1
53
134
208
−1
11
1
Notes for the numerical results based on1000 replications in Table 1:
(1) Bias= the empirical bias × 103;
(2) SE= the empirical standard deviation × 103;
(3) MSE= the empirical mean square error × 104.The outcomes show that the proposed estimator works very well. Through the numerical results in Table1, we derive the following conclusions:(1)
m~(t) outperforms m^(t) obviously in terms of Monte Carlo bias and mean square error for every values of t when n=200 or n=100;(2)
the two estimators are very close to each other in terms of the Monte Carlo standard deviation when the censoring rate is small;(3)
when the censoring rate is heavy,m~(t) outperforms m^(t) obviously in all terms that we referred to in the simulations.Through the conclusions given above, it is appropriate to explore our proposed method in real applications.



---

*Source: 937397-2014-06-09.xml*