Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2744430

Formats

Article sections

- Abstract
- 1. INTRODUCTION
- 2. INFERENCE UNDER INDEPENDENT CENSORING
- 3. INFERENCE UNDER DEPENDENT CENSORING
- 4. EFFICIENCY AND DOUBLE ROBUSTNESS CONSIDERATIONS
- 5. SIMULATION STUDIES
- 6. APPLICATIONS
- 7. CONCLUDING REMARKS
- References

Authors

Related links

J Am Stat Assoc. Author manuscript; available in PMC 2010 June 1.

Published in final edited form as:

J Am Stat Assoc. 2009 June 1; 104(486): 803–815.

doi: 10.1198/jasa.2009.0130PMCID: PMC2744430

NIHMSID: NIHMS132868

Liuquan Sun, Professor, Institute of Applied Mathematics, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100080, P.R. China (E-mail: nc.ca.tma@qls)

See other articles in PMC that cite the published article.

The mean residual life function is an attractive alternative to the survival function or the hazard function of a survival time in practice. It provides the remaining life expectancy of a subject surviving up to time *t*. In this study, we propose a class of transformed mean residual life models for fitting survival data under right censoring. To estimate the model parameters, we make use of the inverse probability of censoring weighting approach and develop a system of estimating equations. Efficiency and robustness of the estimators are also studied. Both asymptotic and finite sample properties of the proposed estimators are established and the approach is applied to two real-life datasets collected from clinical trials.

For a nonnegative survival time *T* with finite expectation, the mean residual life function (MRLF) at time *t* ≥ 0 is defined as

$$m(t)=E(T-t\mid T>t).$$

(1)

In order for the quantity *m*(*t*) to make sense, at least practically, we shall insist that it is only defined when the condition *T* > *t* has positive probability. Throughout this article, we assume *T* is absolutely continuous with a density *f*(*t*) with respect to the Lebesgue measure. Denoting the corresponding survival and hazard functions as *F*(*t*) and *λ*(*t*), respectively, the following four relationships can be easily verified:

$$\begin{array}{l}m(t)=\frac{{\int}_{t}^{\infty}F(u)du}{F(t)},\\ F(t)=\frac{m(0)}{m(t)}exp\left\{-{\int}_{0}^{t}\frac{du}{m(u)}\right\},\\ f(t)=F(t)\left[\frac{\stackrel{.}{m}(t)+1}{m(t)}\right],\end{array}$$

and

$$\lambda (t)=\left[\frac{\stackrel{.}{m}(t)+1}{m(t)}\right],$$

provided that all denominators are nonzero, where (*t*) denotes the first derivative of *m*(*t*).

An important property of MRLF is that *m*(*t*) + *t* is nondecreasing, which is obvious because *m*(*t*) + *t* = *E*(*T*|*T* > *t*). This is also reflected by the fact that (·) + 1 ≥ 0, as *f*(*t*) and *λ*(*t*) are always nonnegative. For other properties of MRLF that do not concern us here, we refer interested readers to Balkema and de Haan (1974), Hollander and Proschan (1975), Kotz and Shanbhag (1980), Arnold and Zahedi (1988), and the references within.

MRLF is of interest in many fields, such as reliability research, survival analysis and actuarial study, and so forth. For example, it is sometimes more informative to tell a prostate cancer patient how long he can survive or live without disease recurrence, in expectation, given his current situation (which of course includes the fact that he has “survived” or lived without the disease so far). As another example, a customer may be interested in knowing how much longer his or her computer can be used, given that the computer has worked normally for, say, *t* years. Later in this article, two real data examples are introduced and analyzed, which will further illustrate the usage of MRLF.

In the following, we will briefly review the literature on MRLF. Early studies on MRLF focused more on its probabilistic behaviors. See, for example, the book by Cox (1961) for a general introduction, and Swartz (1973) and Kotz and Shanbhag (1980) for some characterization and representation theorems. Balkema and de Haan (1974) in particular studied the limiting behavior of *m*(*t*) when *t* → ∞. More recent works include those of Hall and Wellner (1984) and Oakes and Dasu (1990). The concept of MRLF has also been generalized to multivariate distributions. See Arnold and Zahedi (1988) and a recent work by Bassan, Kochar, and Spizzichino (2002).

For statistical inference, many testing procedures have been established. Hollander and Proschan (1975) developed tests for testing exponentiality (under which the MRLF is constant) versus decreasing MRLF alternatives and exponentiality versus “new better than used” in expectation alternatives. In their work, the authors also discussed interesting distribution classes. Recently, El-Bassiouny and Alwasel (2003) proposed a more efficient test for the first scenario considered by Hollander and Proschan (1975). Dauxois (2003), on the other hand, considered testing exponentiality versus nonexponentiality when right-censored data are present.

Estimation of MRLF in homogeneous cases has not been addressed as extensively as the probabilistic properties or testing procedures. Of course, when no censoring is present, this is trivial, because one can easily compute the empirical MRLF. Actually, Bickel et al. (1993, p. 197) showed that it is the optimal (i.e., most efficient) estimator using semiparametric efficiency theory. For bivariate MRLF under order restriction, Rojo and Ghebremichael (2006) developed a projection estimator, again, when no censoring is present. But the literature on censored data seems limited. This issue is discussed briefly in Section 7.

When data are not homogeneous, the estimation problem becomes even more challenging. To our best knowledge, Oakes and Dasu (1990) proposed the first two-sample proportional mean residual life model

$${m}_{1}(t)=\psi {m}_{0}(t),$$

(2)

where *ψ* > 0 represents the ratio of the two MRLFs. They gave a characterization theorem but did not provide a method for estimating *ψ*. Later, in Maguluri and Zhang (1994), model (2) was extended to incorporate covariates,

$$m(t\mid \mathbf{Z})={m}_{0}(t)exp({\mathit{\beta}}^{\prime}\mathbf{Z}),$$

(3)

where *m*(*t*|**Z**) is the MRLF of the survival time *T* corresponding to a *p* × 1 vector covariate **Z**, *m*_{0}(*t*) is some unknown baseline MRLF, and ** β** is an unknown

A new class of additive mean residual life model,

$$m(t\mid \mathbf{Z})={m}_{0}(t)+{\mathit{\beta}}^{\prime}\mathbf{Z},$$

(4)

in which each component has a similar interpretation as their counterparts in model (3), was also proposed by Chen and Cheng (2006) and Chen (2007). In these articles, the authors discussed various estimation methodologies with or without right censoring.

In this article, we propose a more general family of regression models for the MRLF that take the form

$$m(t\mid \mathbf{Z})=g\{{m}_{0}(t)+{\mathit{\beta}}^{\prime}\mathbf{Z}\},$$

(5)

where *g*: ^{+} is a prespecified “transformation” or “link” function, *m*_{0}(*t*) is some unknown function (note that unless *g* is the identity link, *m*_{0}(*t*) is not necessarily the baseline MRLF, whereas *g*{*m*_{0}(*t*)} is), and ** β** and

First, we see that (5) defines a very rich family of models through the link function *g*, which allows the models to be useful in applications such as treatment assessment, resource planning, or short-term forecasting. In particular, it encompasses the Box-Cox transformation, in which *g* is given by *g*(*x*) = [(*x* + 1)* ^{ρ}* − 1]/

Second, we should be careful about the choice of an appropriate link function *g*. Practically, the choice may be based on prior data or the desiring interpretation of the regression parameters. Notice that this choice will also affect the interpretation of the function *m*_{0}(*t*). Theoretically, however, there are some restrictions for choosing an appropriate *g*. In the following, we list the requirements: (A) *g* is twice continuously differentiable, (B) *g* is strictly increasing, and (C) *g*{*m*_{0}(*t*) + ** β**′

Third, as pointed out by a referee, model (5) may prove useful in model selection. Notice that the Box-Cox transformation includes the additive model by letting *ρ* = 1. Thus, if we insist that *g* lies in the family of Box-Cox transformation but allow *ρ* to be estimated, this would provide a way in model selection between the additive mean residual life model and other models contained by the Box-Cox family. We will give a brief discussion in Section 7 on this research topic.

One important motivation behind model (5), in application, is certainly its flexibility in modeling covariate effects. We now look at a real data example. A recent study conducted at the Memorial Sloan-Kettering Cancer Center aims at investigating the potential effect of the radiation therapy dose level on the relapse time of the prostate cancer-specific antigen (PSA). Here, the relapse time, which is our survival time of interest, is defined from the end of the radiation therapy treatment to the time when the PSA level of the patient has been observed to cross a certain threshold value. Of course, this relapse time is sometimes censored by random withdrawal or other reasons. A recent study (Zelefsky et al. 2008) has shown that the age of the patient may also play an important role in predicting the PSA relapse time. Therefore, we need to incorporate age as another covariate in the mean residual life regression model. However, it may not be reasonable to postulate that the effect of age on the MRLF of PSA relapse time is additive or multiplicative, because it would imply that a 5-year increase of age from 60 to 65 would have the same effect as a 5-year increase of age, say, from 80 to 85 (either additively or multiplicatively). As a matter of fact, it is more reasonable clinically to assume that the effect of age phases out gradually as it increases. In other words, we may take *g* to be the Box-Cox transformation with *ρ* as a fraction number like 0.5. More detailed discussion and data analysis results are reported in Section 6.

The rest of this article is organized as follows. In the next section, we will first discuss the situation in which the censoring time is independent of *T* and **Z**, and a general inference procedure based on estimating functions is proposed. The procedure can be easily implemented numerically and the asymptotic properties of the proposed estimators of regression parameters are established. Section 3 extends the method to the situation in which the censoring time may depend on **Z** through the proportional hazards model. Issues regarding efficiency and double robustness are discussed in Section 4. Section 5 reports some results from simulation studies conducted for evaluating the proposed methods. In Section 6, we apply the methodology to two real-life datasets arising from clinical trials, and some concluding remarks are given in Section 7. Regularity conditions and theoretical proofs are collected in the Appendix and a supplemental technical report, respectively.

Let *C* be the potential censoring time, and assume that the support of *C* is longer than that of the survival time *T* to ensure that the MRLF is estimable. In the technical arguments, this assumption also saves us from lengthy technical discussion of the tail behavior of the limiting distributions. In this section, we consider the situation when *C* is independent of both *T* and **Z**. Let {*T _{i}*,

Let *G*(*t*) be the survival function of *C* and *Ĝ*(*t*) the Kaplan-Meier estimator of *G*(*t*) based on {*X _{i}*, 1 − Δ

$$\begin{array}{l}{M}_{i}(t)=\frac{{\mathrm{\Delta}}_{i}I({X}_{i}>t)}{G({X}_{i})}[({X}_{i}-t)-g\{{m}_{0}(t)+{\mathit{\beta}}_{0}^{\prime}{\mathbf{Z}}_{i}\}],\\ i=1,\dots ,n.\end{array}$$

Under model (5), *M _{i}*(

$$\sum _{i=1}^{n}\frac{{\mathrm{\Delta}}_{i}I({X}_{i}>t)}{\widehat{G}({X}_{i})}[({X}_{i}-t)-g\{{m}_{0}(t)+{\mathit{\beta}}^{\prime}{\mathbf{Z}}_{i}\}]=0,0\le t\le \tau ,$$

(6)

where 0 < τ = inf {*t: Pr*(*T* ≥ *t*) = 0} < ∞. Note that, without censoring, (6) reduces to
${\sum}_{i=1}^{n}I({T}_{i}>t)[({T}_{i}-t)-g\{{m}_{0}(t)+{\mathit{\beta}}^{\prime}{\mathbf{Z}}_{i}\}]=0$, which is similar to the efficient nonparametric estimator of the MRLF discussed in Bickel et al. (1993). The inverse probability of censoring weighting (IPCW) technique has been extensively discussed in the literature. See, for example, Robins and Rotnitzky (1992); Robins, Rotnitzky, and Zhao (1994); and van der Laan and Robins (2003).

Denote this estimator by _{a}_{0}(*t*; ** β**). To estimate

$$\begin{array}{l}{\mathbf{U}}_{a}(\mathit{\beta})=\sum _{i=1}^{n}{\int}_{0}^{\tau}\frac{{\mathrm{\Delta}}_{i}I({X}_{i}>t){\mathbf{Z}}_{i}}{\widehat{G}({X}_{i})}[({X}_{i}-t)-g\{{\widehat{m}}_{a0}(t;\mathit{\beta})\\ +{\mathit{\beta}}^{\prime}{\mathbf{Z}}_{i}\}]dH(t)=\mathbf{0},\end{array}$$

(7)

where *H*(*t*) is an increasing weight function on [0, *τ*]. Let * _{a}* denote the solution to

$$\begin{array}{l}{\widehat{M}}_{i}^{c}(t)={N}_{i}^{c}(t)-{\int}_{0}^{t}I({X}_{i}\ge u)d{\widehat{\mathrm{\Lambda}}}^{c}(u),\\ {\widehat{M}}_{i}(t)=\frac{{\mathrm{\Delta}}_{i}I({X}_{i}>t)}{\widehat{G}({X}_{i})}\left[({X}_{i}-t)-g\{{\widehat{m}}_{a0}(t)+{{\widehat{\mathit{\beta}}}^{\prime}}_{a}{\mathbf{Z}}_{i}\}\right],\\ {\widehat{\mathbf{Q}}}_{a}(t)={n}^{-1}\sum _{i=1}^{n}I({X}_{i}\ge t){\int}_{0}^{\tau}{\widehat{M}}_{i}(u)\left\{{\mathbf{Z}}_{i}-{\overline{\mathbf{Z}}}_{a}(u;{\widehat{\mathit{\beta}}}_{a})\right\}dH(u)\end{array}$$

and

$${\overline{\mathbf{Z}}}_{a}(t;\mathit{\beta})=\frac{{\sum}_{i=1}^{n}{\mathrm{\Delta}}_{i}\widehat{G}{({X}_{i})}^{-1}I({X}_{i}>t)\stackrel{.}{g}\{{\widehat{m}}_{a0}(t;\mathit{\beta})+{\mathit{\beta}}^{\prime}{\mathbf{Z}}_{i}\}{\mathbf{Z}}_{i}}{{\sum}_{i=1}^{n}{\mathrm{\Delta}}_{i}\widehat{G}{({X}_{i})}^{-1}I({X}_{i}>t)\stackrel{.}{g}\{{\widehat{m}}_{a0}(t;\mathit{\beta})+{\mathit{\beta}}^{\prime}{\mathbf{Z}}_{i}\}},$$

where *ġ*(*x*) = *dg*(*x*)/*dx*. We summarize the asymptotic properties of * _{a}* and

Under the regularity conditions (C1) through (C4) stated in the Appendix, we have the following:

and_{a}_{a}_{0}(*t*) always exist, and are unique and consistent.*n*^{1/2}(−_{a}*β*_{0}) is asymptotically normal with mean zero and a variance–covariance matrix that can be consistently estimated by**Â**^{−1}**Σ**_{a}**Â**^{−1}, where ${\mathrm{\sum}}_{a}={n}^{-1}{\sum}_{i=1}^{n}{\widehat{\mathit{\xi}}}_{i}^{\otimes 2}$,$$\begin{array}{l}{\widehat{\mathit{\xi}}}_{i}={\int}_{0}^{\tau}{\widehat{M}}_{i}(t)\left\{{\mathbf{Z}}_{i}-{\overline{\mathbf{Z}}}_{a}(t;{\widehat{\mathit{\beta}}}_{a})\right\}dH(t)+{\int}_{0}^{\tau}\frac{{\widehat{\mathbf{Q}}}_{a}(t)}{\widehat{\pi}(t)}d{\widehat{M}}_{i}^{c}(t),\\ \widehat{\mathbf{A}}={n}^{-1}\sum _{i=1}^{n}{\int}_{0}^{\tau}\frac{{\mathrm{\Delta}}_{i}I({X}_{i}>t)}{\widehat{G}({X}_{i})}{\{{\mathbf{Z}}_{i}-{\overline{\mathbf{Z}}}_{a}(t;{\widehat{\mathit{\beta}}}_{a})\}}^{\otimes 2}\stackrel{.}{g}\{{\widehat{m}}_{a0}(t)\\ +{\widehat{\mathit{\beta}}}_{a}^{\prime}{\mathbf{Z}}_{i}\}dH(t)\end{array}$$and**v**^{2}=**vv**′ for a column vector**v**.*n*^{1/2}{_{a}_{0}(*t*) −*m*_{0}(*t*)}(0 ≤*t*≤ τ) converges weakly to a zero-mean Gaussian process with a covariance function that can be consistently estimated by ${\widehat{\mathrm{\Gamma}}}_{a}(s,t)={n}^{-1}{\sum}_{i=1}^{n}{\widehat{\varphi}}_{i}(s){\widehat{\varphi}}_{i}(t)$at (*s, t*), where$$\begin{array}{l}{\widehat{\varphi}}_{i}(t)=\frac{1}{{\widehat{\mathrm{\Phi}}}_{a}(t)}\left[{\widehat{M}}_{i}(t)+{\int}_{0}^{\tau}\frac{{\widehat{R}}_{a}(t,u)}{\widehat{\pi}(u)}d{\widehat{M}}_{i}^{c}(u)\right]-{\overline{\mathbf{Z}}}_{a}(t;{\widehat{\mathit{\beta}}}_{a}{)}^{\prime}{\widehat{\mathbf{A}}}^{-1}{\widehat{\mathit{\xi}}}_{i},\\ {\widehat{\mathrm{\Phi}}}_{a}(t)={n}^{-1}\sum _{i=1}^{n}\frac{{\mathrm{\Delta}}_{i}I({X}_{i}>t)}{\widehat{G}({X}_{i})}\stackrel{.}{g}\{{\widehat{m}}_{a0}(t)+{{\widehat{\mathit{\beta}}}^{\prime}}_{a}{\mathbf{Z}}_{i}\},\end{array}$$and ${\widehat{R}}_{a}(t,u)={n}^{-1}{\sum}_{i=1}^{n}{\widehat{M}}_{i}(t)I({X}_{i}\ge u)$.

The asymptotic normality for _{a}_{0}(*t*), together with the consistent variance estimator * _{a}*, enables us to construct pointwise confidence intervals for

We now consider the situation in which the censoring time *C* may depend on the covariate **Z**. However, we assume that given **Z**, *T* is independent of *C*. This actually corresponds to the scenario called “coarsening at random” in the literature of missing data. It has been pointed out (Robins and Ritov 1997) that for such problems it is wise to postulate a model for the censoring mechanism. Of course, if **Z** takes only finitely many values and the sample size *n* is large enough, we may be able to calculate the Kaplan-Meier estimates of the censoring distribution for each value of **Z**. However, this is not feasible when **Z** (or at least one component of **Z**) is continuous or the sample size is moderate. To this end, we assume that the hazard function of *C* given **Z** has the form

$${\lambda}_{c}(t\mid \mathbf{Z})={\lambda}_{0}(t)exp\{{\gamma}^{\prime}\mathbf{Z}\},$$

(8)

where *λ*_{0}(*t*) is an unspecified baseline hazard function and ** γ** is a vector of unknown regression parameters. Similar as before, we denote the true value of

Note that the model from Section 2 is the model here with *γ*_{0} = **0**. Of course, *γ*_{0}, as well as *λ*_{0}(*t*), needs to be estimated. A natural estimator of *γ*_{0} under model (8) is given by the convenient maximum partial likelihood estimator defined as the solution to

$${\mathbf{U}}_{r}(\mathit{\gamma})=\sum _{i=1}^{n}{\int}_{0}^{\tau}\{{\mathbf{Z}}_{i}-{\overline{\mathbf{Z}}}_{r}(t;\mathit{\gamma})\}d{N}_{i}^{c}(t)=\mathbf{0},$$

(9)

where
${N}_{i}^{c}(t)=I({X}_{i}\le t,{\mathrm{\Delta}}_{i}=0)$, * _{r}*(

$${\widehat{\mathrm{\Lambda}}}_{0}(t)=\sum _{i=1}^{n}{\int}_{0}^{t}\frac{d{N}_{i}^{c}(u)}{{\sum}_{j=1}^{n}I({X}_{j}\ge u)exp\{{\widehat{\mathit{\gamma}}}^{\prime}{\mathbf{Z}}_{j}\}}.$$

Also, let *G _{i}*(

$$\begin{array}{l}\sum _{i=1}^{n}\frac{{\mathrm{\Delta}}_{i}I({X}_{i}>t)}{{G}_{i}({X}_{i};\widehat{\mathit{\gamma}},{\widehat{\mathrm{\Lambda}}}_{0})}[({X}_{i}-t)-g\{{m}_{0}(t)+{\mathit{\beta}}^{\prime}{\mathbf{Z}}_{i}\}]=0,\\ 0\le t\le \tau .\end{array}$$

(10)

Denote this estimator by _{b}_{0}(*t*; ** β**) for a given

$$\begin{array}{l}{\mathbf{U}}_{b}(\mathit{\beta})\\ =\sum _{i=1}^{n}{\int}_{0}^{\tau}\frac{{\mathrm{\Delta}}_{i}I({X}_{i}>t){\mathbf{Z}}_{i}}{{G}_{i}({X}_{i};\widehat{\mathit{\gamma}},{\widehat{\mathrm{\Lambda}}}_{0})}[({X}_{i}-t)-g\{{\widehat{m}}_{b0}(t;\mathit{\beta})+{\mathit{\beta}}^{\prime}{\mathbf{Z}}_{i}\}]\\ \times dH(t)=\mathbf{0}.\end{array}$$

(11)

Let * _{b}* denote the solution to

$$\begin{array}{l}{\stackrel{\sim}{M}}_{i}^{c}(t)={N}_{i}^{c}(t)-{\int}_{0}^{t}I({X}_{i}\ge u)exp\{{\widehat{\mathit{\gamma}}}^{\prime}{\mathbf{Z}}_{i}\}d{\widehat{\mathrm{\Lambda}}}_{0}(u),\\ {\widehat{M}}_{i}^{\ast}(t)=\frac{{\mathrm{\Delta}}_{i}I({X}_{i}>t)}{{G}_{i}({X}_{i};\widehat{\mathit{\gamma}},{\widehat{\mathrm{\Lambda}}}_{0})}\left[({X}_{i}-t)-g\{{\widehat{m}}_{b0}(t)+{\widehat{\mathit{\beta}}}_{b}^{\prime}{\mathbf{Z}}_{i}\}\right],\\ {\widehat{\mathbf{Q}}}_{b}(t)={n}^{-1}\sum _{i=1}^{n}I({X}_{i}\ge t)exp\{{\widehat{\mathit{\gamma}}}^{\prime}{\mathbf{Z}}_{i}\}\\ \times {\int}_{0}^{\tau}{\widehat{M}}_{i}^{\ast}(u)\left\{{\mathbf{Z}}_{i}-{\overline{\mathbf{Z}}}_{b}(u;{\widehat{\mathit{\beta}}}_{b})\right\}dH(u),\\ \widehat{\mathbf{P}}={n}^{-1}\sum _{i=1}^{n}{\widehat{\mathrm{\Lambda}}}_{0}({X}_{i})exp\{{\widehat{\mathit{\gamma}}}^{\prime}{\mathbf{Z}}_{i}\}\\ \times {\int}_{0}^{\tau}{\widehat{M}}_{i}^{\ast}(u)\left\{{\mathbf{Z}}_{i}-{\overline{\mathbf{Z}}}_{b}(u;{\widehat{\mathit{\beta}}}_{b})\right\}dH(u){{\mathbf{Z}}^{\prime}}_{i},\\ \frac{{\overline{\mathbf{Z}}}_{b}(t;\mathit{\beta})={\sum}_{i=1}^{n}{\mathrm{\Delta}}_{i}{G}_{i}{({X}_{i};\widehat{\mathit{\gamma}},{\widehat{\mathrm{\Lambda}}}_{0})}^{-1}I({X}_{i}>t)\stackrel{.}{g}\{{\widehat{m}}_{b0}(t;\mathit{\beta})+{\mathit{\beta}}^{\prime}{\mathbf{Z}}_{i}\}{\mathbf{Z}}_{i}}{{\sum}_{i=1}^{n}{\mathrm{\Delta}}_{i}{G}_{i}{({X}_{i};\widehat{\mathit{\gamma}},{\widehat{\mathrm{\Lambda}}}_{0})}^{-1}I({X}_{i}>t)\stackrel{.}{g}\{{\widehat{m}}_{b0}(t;\mathit{\beta})+{\mathit{\beta}}^{\prime}{\mathbf{Z}}_{i}\}}\end{array}$$

We summarize the asymptotic properties of * _{b}* and

Under the regularity conditions (C1) through (C3) and (C5) stated in the Appendix, we have

and_{b}_{b}_{0}(*t*) always exist, and are unique and consistent.*n*^{1/2}(−_{b}*β*_{0}) is asymptotically normal with mean zero and a variance–covariance matrix that can be consistently estimated by^{−1}_{b}^{−1}, where ${\widehat{\mathbf{\sum}}}_{b}={n}^{-1}{\sum}_{i=1}^{n}{\widehat{\mathit{\eta}}}_{i}^{\otimes 2}$,$$\begin{array}{l}{\widehat{\mathit{\eta}}}_{i}={\int}_{0}^{\tau}{\widehat{M}}_{i}^{\ast}(t)\left\{{\mathbf{Z}}_{i}-{\overline{\mathbf{Z}}}_{b}(t;{\widehat{\mathit{\beta}}}_{b})\right\}dH(t)+{\int}_{0}^{\tau}\frac{{\widehat{\mathbf{Q}}}_{b}(t)}{{S}^{(0)}(t;\widehat{\mathit{\gamma}})}d{\stackrel{\sim}{M}}_{i}^{c}(t)\\ +\left[\widehat{\mathbf{P}}-{\int}_{0}^{\tau}{\widehat{\mathbf{Q}}}_{b}(t){\overline{\mathbf{Z}}}_{r}(t;\widehat{\mathit{\gamma}}{)}^{\prime}d{\widehat{\mathrm{\Lambda}}}_{0}(t)\right]{\widehat{\mathbf{D}}}^{-1}\\ \times {\int}_{0}^{\tau}\{{\mathbf{Z}}_{i}-{\overline{\mathbf{Z}}}_{r}(t;\widehat{\mathit{\gamma}})\}d{\stackrel{\sim}{M}}_{i}^{c}(t),\\ \widehat{\mathbf{B}}={n}^{-1}\sum _{i=1}^{n}{\int}_{0}^{\tau}\frac{{\mathrm{\Delta}}_{i}I({X}_{i}>t)}{{G}_{i}({X}_{i};\widehat{\mathit{\gamma}},{\widehat{\mathrm{\Lambda}}}_{0})}{\{{\mathbf{Z}}_{i}-{\overline{\mathbf{Z}}}_{b}(t;{\widehat{\mathit{\beta}}}_{b})\}}^{\otimes 2}\stackrel{.}{g}\{{\widehat{m}}_{b0}(t)\\ +{\widehat{\mathit{\beta}}}_{b}^{\prime}{\mathbf{Z}}_{i}\}dH(t),\end{array}$$and = −*n*^{−1}**U**(_{r}**)/***γ*′.*n*^{1/2}{_{b}_{0}(*t*) −*m*_{0}(*t*)} (0≤*t*≤*τ*) converges weakly to a zero-mean Gaussian process with a covariance function at (*s, t*) that can be consistently estimated by ${\widehat{\mathrm{\Gamma}}}_{b}(s,t)={n}^{-1}{\sum}_{i=1}^{n}{\widehat{\psi}}_{i}(s){\widehat{\psi}}_{i}(t)$, where$$\begin{array}{l}{\widehat{\psi}}_{i}(t)=\frac{1}{{\widehat{\mathrm{\Phi}}}_{b}(t)}\left[{\widehat{M}}_{i}^{\ast}(t)+{\int}_{0}^{\tau}\frac{{\widehat{R}}_{b}(t,u)}{{S}^{(0)}(u;\widehat{\mathit{\gamma}})}d{\stackrel{\sim}{M}}_{i}^{c}(u)\right]\\ -{\overline{\mathbf{Z}}}_{b}(t;{\widehat{\mathit{\beta}}}_{b}{)}^{\prime}{\widehat{\mathbf{B}}}^{-1}{\widehat{\mathit{\eta}}}_{i}+\widehat{\mathit{\psi}}(t){D}^{-1}{\int}_{0}^{\tau}\{{\mathbf{Z}}_{i}-{\overline{\mathbf{Z}}}_{r}(u;\widehat{\mathit{\gamma}})\}d{\stackrel{\sim}{M}}_{i}^{c}(u),\\ {\widehat{\mathrm{\Phi}}}_{b}(t)={n}^{-1}\sum _{i=1}^{n}\frac{{\mathrm{\Delta}}_{i}I({X}_{i}>t)}{{G}_{i}({X}_{i};\widehat{\mathit{\gamma}},{\widehat{\mathrm{\Lambda}}}_{0})}\stackrel{.}{g}\{{\widehat{m}}_{b0}(t)+{\widehat{\mathit{\beta}}}_{b}^{\prime}{\mathbf{Z}}_{i}\},\\ \widehat{\mathit{\psi}}(t)={n}^{-1}\sum _{i=1}^{n}{\widehat{M}}_{i}^{\ast}(t){\widehat{\mathrm{\Lambda}}}_{0}({X}_{i})exp\{{\widehat{\mathit{\gamma}}}^{\prime}{\mathbf{Z}}_{i}\}\\ \times {{\mathbf{Z}}^{\prime}}_{i}-{\int}_{0}^{t}{\widehat{R}}_{b}(t,u){\overline{\mathbf{Z}}}_{r}(u;\widehat{\mathit{\gamma}}{)}^{\prime}d{\widehat{\mathrm{\Lambda}}}_{0}(u),\\ {\widehat{R}}_{b}(t,u)={n}^{-1}\sum _{i=1}^{n}{\widehat{M}}_{i}^{\ast}(t)exp\{{\widehat{\mathit{\gamma}}}^{\prime}{\mathbf{Z}}_{i}\}I({X}_{i}\ge u).\end{array}$$

To construct a simultaneous confidence band for *m*_{0}(*t*) over [*t*_{1}, *t*_{2}], we may again use the resampling method mentioned in Section 2. Namely, we approximate the distribution of the process *n*^{1/2}{_{b}_{0}(*t*) − *m*_{0}(*t*)} by that of the zero-mean Gaussian process *Ŵ _{b}*(

Because estimating functions (7) and (11) are given in a somewhat ad hoc fashion using the IPCW method, a question naturally arises regarding how efficient and robust the resulting inference procedures are. Here we provide insights into these problems by using the semiparametric efficiency theory developed by Bickel et al. (1993); Robins, Rotnitzky, and Zhao (1994); and van der Laan and Robins (2003). For simplicity, we will only study the efficiency and double robustness in the case of independent censoring. A discussion on the case of dependent censoring, although similar to those given later conceptually, is much messier algebraically and is thus omitted.

The fundamental idea of semiparametric efficiency is conveyed through the so-called “tangent spaces” and orthogonal projections in Hilbert spaces. Normally, one makes enough regularity assumptions in order for the problem to be meaningful practically and solvable theoretically (Newey, 1990). For instance, it is usually insisted that we deal only with regular and asymptotically linear estimators. This excludes superefficient estimators and allows us to identify influence functions. Notice that this is exactly the situation we have here for * _{a}* (see Theorem 1 and its proof). Under these assumptions, the tangent space is taken to be the closed linear span of the scores of the parametric submodels of the semiparametric model of interest. One may then orthogonally project an arbitrary influence function onto this tangent space and obtain the efficient influence function. Alternatively, if the parametric part and the nonparametric part of the semiparametric model can be explicitly separated, it might be possible to express the tangent space by the direct sum of two subspaces: one corresponding to the parameter of interest and the other corresponding to the nuisance parameter. The latter is sometimes referred to as the “nuisance tangent space.” One may then take the orthogonal projection of the score of the parametric part onto the nuisance tangent space and obtain the efficient score function.

These approaches, although theoretically elegant, may not be computationally feasible. The difficulty usually lies in two aspects: (1) finding the tangent space (or the nuisance tangent space) and (2) calculating the orthogonal projection. The first difficulty is mostly the result of the complicated structure of the semiparametric model, and thus the scores. The second difficulty may be the result of several reasons. For example, the orthogonal projection might be a conditional expectation that requires extra modeling and awkward numerical integration. When data are coarsened or censored and the coarsening or censoring mechanism needs to be modeled, these difficulties multiply. For our problem, unfortunately, both difficulties exist. Therefore, we take the following approach. First, because * _{a}* has already been shown to be regular and asymptotically linear, we will take advantage of its estimating function

First, notice that we have already derived the influence function for * _{a}* in Theorem 1 and its proof. Namely,

$$\begin{array}{l}{n}^{1/2}({\widehat{\mathit{\beta}}}_{a}-{\mathit{\beta}}_{0})={\mathbf{A}}^{-1}{n}^{-1/2}{\mathbf{U}}_{a}({\mathit{\beta}}_{0})+{o}_{p}(1)\\ ={\mathbf{A}}^{-1}{n}^{-1/2}\sum _{i=1}^{n}{\mathit{\xi}}_{i}+{o}_{p}(1).\end{array}$$

In other words, the influence function is **A**^{−1}*ξ** _{i}*. For a definition of

$$\begin{array}{l}{\mathbf{D}}_{i}={\int}_{0}^{\tau}I({T}_{i}>t)[({T}_{i}-t)-g\{{m}_{0}(t)+{\mathit{\beta}}_{0}^{\prime}{\mathbf{Z}}_{i}\}]\{{\mathbf{Z}}_{i}-{\overline{\mathbf{z}}}_{a}(t)\}\\ \times d\stackrel{\sim}{H}(t),\phantom{\rule{0.38889em}{0ex}}i=1,\dots ,n,\end{array}$$

and **D** accordingly, where * _{a}*(

$${M}_{i}^{c}(t)={N}_{i}^{c}(t)-{\int}_{0}^{t}I({X}_{i}\ge u)d{\mathrm{\Lambda}}^{c}(u),\phantom{\rule{0.38889em}{0ex}}i=1,\dots ,n.$$

Then
${M}_{i}^{c}(t)$are martingales with respect to the *σ*-filtration

$$\begin{array}{l}\sigma \{\text{I}({X}_{i}\ge u),I({X}_{i}\le u,{\mathrm{\Delta}}_{i}=0),0\le u\le t;I({T}_{i}\le \upsilon ),\\ 0\le \upsilon \le \tau ;{\mathbf{Z}}_{i},i=1,\dots ,n\}.\end{array}$$

Similar settings can be found in Zhao and Tsiatis (1999) and Bang and Tsiatis (2000, 2002). Consider the martingale representation of the Kaplan-Meier estimator (Fleming and Harrington, 1991, p. 97)

$$\frac{\widehat{G}(t)-G(t)}{G(t)}=-{\int}_{0}^{t}\frac{\widehat{G}(u-)}{G(u)}\frac{{\sum}_{i=1}^{n}d{M}_{i}^{c}(u)}{n\widehat{\pi}(u)}$$

(12)

and the result from Robins and Rotnitzky (1992) that

$$\frac{{\mathrm{\Delta}}_{i}}{G({X}_{i})}=1-{\int}_{0}^{\tau}\frac{d{M}_{i}^{c}(t)}{G(t)}.$$

(13)

It then follows from (12), (13), and some algebraic manipulation that

$$\begin{array}{l}{\mathbf{U}}_{a}({\beta}_{0})=\sum _{i=1}^{n}\frac{{\mathrm{\Delta}}_{i}{\mathbf{D}}_{i}}{G({X}_{i})}-\sum _{i=1}^{n}{\mathrm{\Delta}}_{i}{\mathbf{D}}_{i}\frac{\widehat{G}({X}_{i})-G({X}_{i})}{G({X}_{i})\widehat{G}({X}_{i})}+{o}_{p}({n}^{1/2})\\ =\sum _{i=1}^{n}\left[{\mathbf{D}}_{i}-{\int}_{0}^{\tau}\{{\mathbf{D}}_{i}-\mathbf{K}(t;\mathbf{D})\}\frac{d{M}_{i}^{c}(t)}{G(t)}\right]+{o}_{p}({n}^{1/2}),\end{array}$$

where *K*(*t; W*) = *E*{*WI*(*T* ≥ *t*)}/*F*(*t*) for any random variable or functional *W*, and *F*(*t*) = *P*{*T* ≥ *t*} is the (marginal) survival function of *T*. This means the *i*th influence function of the regular and asymptotically linear estimator * _{a}* is

$${\mathbf{A}}^{-1}\left[{\mathbf{D}}_{i}-{\int}_{0}^{\tau}\{{\mathbf{D}}_{i}-\mathbf{K}(t;\mathbf{D})\}\frac{d{M}_{i}^{c}(t)}{G(t)}\right].$$

(14)

The alternative expression (14) for **U*** _{a}*(

We now seek an AIPCW estimator that is more efficient than * _{a}*. Put another way, it should have a smaller asymptotic variance than (14). Here, by convention, we say that a variance matrix

$$\begin{array}{l}{\mathbf{D}}_{i}-{\int}_{0}^{\tau}\{{\mathbf{D}}_{i}-\mathbf{K}(t;\mathbf{D})\}\frac{d{M}_{i}^{c}(t)}{G(t)}\\ +{\int}_{0}^{\tau}\{\mathbf{e}(t;{\mathbf{Z}}_{i})-\mathbf{K}(t;\mathbf{e})\}\frac{d{M}_{i}^{c}(t)}{G(t)},\end{array}$$

(15)

where **e**(*t*; **Z*** _{i}*) is an arbitrary

$${\mathbf{V}}_{1}=E\left[{\int}_{0}^{\tau}{\mathbf{D}}_{i}\{\mathbf{e}(t;{\mathbf{Z}}_{i})-\mathbf{K}(t;\mathbf{e}){\}}^{\prime}I({T}_{i}\ge t){\scriptstyle \frac{d{\mathrm{\Lambda}}^{c}(t)}{G(t)}}\right].$$

and

$${\mathbf{V}}_{2}=E\left[{\int}_{0}^{\tau}{\{\mathbf{e}(t;{\mathbf{Z}}_{i})-\mathbf{K}(t;\mathbf{e})\}}^{\otimes 2}I({X}_{i}\ge t){\scriptstyle \frac{d{\mathrm{\Lambda}}^{c}(t)}{G{(t)}^{2}}}\right].$$

Note that **V**_{1} is actually the covariance of the second and the third terms in (15) whereas **V**_{2} is the variance of the third term using the stochastic integration theory. With this modification, the “improved” function (15) now becomes

$$\begin{array}{l}{\mathbf{D}}_{i}-{\int}_{0}^{\tau}\{{\mathbf{D}}_{i}-\mathbf{K}(t;\mathbf{D})\}\frac{d{M}_{i}^{c}(t)}{G(t)}\\ +\mathrm{\Upsilon}{\int}_{0}^{\tau}\{\mathbf{e}(t;{\mathbf{Z}}_{i})-\mathbf{K}(t;\mathbf{e})\}\frac{d{M}_{i}^{c}(t)}{G(t)},\end{array}$$

(16)

with a variance that is

$$\begin{array}{l}var\{{\mathbf{D}}_{i}\}+E\left[{\int}_{0}^{\tau}{\{{\mathbf{D}}_{i}-\mathbf{K}(t;\mathbf{D})\}}^{\otimes 2}I({X}_{i}\ge t)\frac{d{\mathrm{\Lambda}}^{c}(t)}{G{(t)}^{2}}\right]\\ -{\mathbf{V}}_{1}{\mathbf{V}}_{2}^{-1}{{\mathbf{V}}^{\prime}}_{1},\end{array}$$

(17)

which is always smaller than that of the term in the bracket of (14). This implies that if we can construct an estimating equation based on (16), then the resulting estimator should have a smaller variance than * _{a}* as a result of (17) and the fact that the augmented term in (16) is independent of

$${\stackrel{\sim}{\mathbf{D}}}_{i}(\mathit{\beta})={\int}_{0}^{\tau}I({T}_{i}>t){\mathbf{Z}}_{i}[({T}_{i}-t)-g\{{\widehat{m}}_{a0}(t;\mathit{\beta})+{\mathit{\beta}}^{\prime}{\mathbf{Z}}_{i}\}]dH(t).$$

A class of AIPCW estimating equations according to (16) is given by

$${\stackrel{\sim}{\mathbf{U}}}_{a}(\mathit{\beta})=\sum _{i=1}^{n}\{\frac{{\mathrm{\Delta}}_{i}{\stackrel{\sim}{\mathbf{D}}}_{i}(\mathit{\beta})}{\widehat{G}({X}_{i})}+\widehat{\mathrm{\Upsilon}}({\widehat{\mathit{\beta}}}_{a}){\int}_{0}^{\tau}\{\mathbf{e}(t;{\mathbf{Z}}_{i})-\widehat{\mathbf{K}}(t;\mathbf{e})\}\frac{d{N}_{i}^{c}(t)}{\widehat{G}(t)}\},$$

(18)

where $\widehat{\mathbf{K}}(t;\mathbf{e})={\sum}_{i=1}^{n}\mathbf{e}(t;{\mathbf{Z}}_{i})I({X}_{i}\ge t)/{\sum}_{i=1}^{n}I({X}_{i}\ge t),\widehat{\mathrm{\Upsilon}}({\widehat{\mathit{\beta}}}_{a})={\widehat{\mathbf{V}}}_{1}({\widehat{\mathit{\beta}}}_{a}){\widehat{\mathbf{V}}}_{2}^{-1}$,

$$\begin{array}{l}{\widehat{\mathbf{V}}}_{1}({\widehat{\mathit{\beta}}}_{a})={n}^{-1}\sum _{i=1}^{n}{\int}_{0}^{\tau}\frac{{\mathrm{\Delta}}_{i}{\stackrel{\sim}{\mathbf{D}}}_{i}({\widehat{\mathit{\beta}}}_{a})}{\widehat{G}({X}_{i})}\{\mathbf{e}(t;{\mathbf{Z}}_{i})-\widehat{\mathbf{K}}(t;\mathbf{e}){\}}^{\prime}I({X}_{i}\ge t)\\ \times \frac{d{\widehat{\mathrm{\Lambda}}}^{c}(t)}{\widehat{G}(t)},\end{array}$$

and

$${\widehat{\mathbf{V}}}_{2}={n}^{-1}\sum _{i=1}^{n}{\int}_{0}^{\tau}{\{\mathbf{e}(t;{\mathbf{Z}}_{i})-\widehat{\mathbf{K}}(t;\mathbf{e})\}}^{\otimes 2}I({X}_{i}\ge t)\frac{d{\widehat{\mathrm{\Lambda}}}^{c}(t)}{\widehat{G}{(t)}^{2}}.$$

Let * _{a}* be the improved estimator—in other words, the solution to

$$\begin{array}{l}{\mathbf{\sum}^{\sim}}_{a}={n}^{-1}\sum _{i=1}^{n}\frac{{\mathrm{\Delta}}_{i}{\stackrel{\sim}{\mathbf{D}}}_{i}{({\stackrel{\sim}{\mathit{\beta}}}_{a})}^{\otimes 2}}{\widehat{G}({X}_{i})}+{n}^{-1}\sum _{i=1}^{n}{\int}_{0}^{\tau}\{\stackrel{\sim}{\mathbf{K}}(t;{\stackrel{\sim}{\mathbf{D}}}^{\otimes 2})\\ -\stackrel{\sim}{\mathbf{K}}{(t;\stackrel{\sim}{\mathbf{D}})}^{\otimes 2}\}\frac{d{N}_{i}^{c}(t)}{\widehat{G}{(t)}^{2}}-{\widehat{\mathbf{V}}}_{1}({\widehat{\mathit{\beta}}}_{a}){\widehat{\mathbf{V}}}_{2}^{-1}{\widehat{\mathbf{V}}}_{1}({\widehat{\mathit{\beta}}}_{a}{)}^{\prime}\mid \end{array}$$

and
$\stackrel{\sim}{\mathbf{K}}(t;\mathbf{W})={\{n\widehat{F}(t)\}}^{-1}{\sum}_{i=1}^{n}{\mathrm{\Delta}}_{i}\mathbf{W}I({X}_{i}\ge t)/\widehat{G}({X}_{i})$for any random vector or matrix **W**.

This class of improved estimators is generated through the functional of **e**(*t*; **Z*** _{i}*). In practice, we can obtain an improved estimator by choosing any form of

The augmented estimating function **Ũ*** _{a}*(

In a censoring data problem, an estimator is double robust (DR) if it remains consistent when either the model for the censoring mechanism or a model for the complete data distribution is correctly specified (see Robins and Rotnitzky (2001) for an existence theorem). For our problem, we wish to provide such “double protection” because the censoring distribution might be wrongly modeled. In general, double robustness may result in certain efficiency lost; however, the AIPCW estimator does not have this shortcoming. As argued in a series of studies (e.g., Scharfstein, Rotnitzky, and Robins 1999; van der Laan and Robins 2003; Rotnitzky and Robins 2005), the optimal AIPCW estimator provides both efficiency gain and double robustness.

In the following, we consider the DR property of the optimal estimator derived in the last subsection. The major mathematical tools can be found in van der Laan and Robins (2003, chaps. 1 and 2) and Tsiatis (2006, section 10.3). Recall we have concluded that the efficient influence function is obtained if we take **e**(*t*; **Z*** _{i}*) =

$$\frac{{\mathrm{\Delta}}_{i}{\mathbf{D}}_{i}}{G({X}_{i})}+{\int}_{0}^{\tau}E[{\mathbf{D}}_{i}\mid {T}_{i}\ge t,{\mathbf{Z}}_{i}]\frac{d{M}_{i}^{c}(t)}{G(t)}.$$

(19)

For notational convenience we will now write **D*** _{i}* as

$$\begin{array}{l}{\mathbf{Q}}_{F}(t,\mathbf{Z}):=-E[\mathbf{D}(T,\mathbf{Z})\mid T\ge t,\mathbf{Z}]\\ =\frac{1}{F(t\mid \mathbf{Z})}{\int}_{t}^{\tau}\mathbf{D}(u,\mathbf{Z})dF(u\mid \mathbf{Z}),\end{array}$$

where *F*(*t*|**Z**) is the survival function of *T* given **Z**. Then (19) becomes

$${\mathbf{Y}}_{F,G}(X,\mathrm{\Delta},\mathbf{Z}):=\frac{\mathrm{\Delta}\mathbf{D}(X,\mathbf{Z})}{G(X)}-{\int}_{0}^{\tau}{\mathbf{Q}}_{F}(t,\mathbf{Z})\frac{d{M}^{c}(t)}{G(t)}.$$

(20)

Equation (20) implies that the efficient influence function of the estimator depends on two distribution functions, *F*(·|**Z**) and *G*(·), which, in practice, will both need to be modeled (in contrast, the simple IPCW estimator only requires modeling *G*(·), which under independent censoring assumption is the nonparametric Kaplan-Meier estimator). As in Rubin and van der Laan (2006), where **Y**_{F}_{,}* _{G}*(

$$E\{{\mathbf{Y}}_{{F}_{1},{G}_{1}}(X,\mathrm{\Delta},\mathbf{Z})\mid \mathbf{Z}\}=\mathbf{0}\phantom{\rule{0.38889em}{0ex}}\text{if}\phantom{\rule{0.16667em}{0ex}}{F}_{1}=F\phantom{\rule{0.38889em}{0ex}}\text{or}\phantom{\rule{0.16667em}{0ex}}{G}_{1}=G.$$

(21)

This means that we only need to model one of them correctly. That is, we get two chances of obtaining a consistent estimator. This double robustness property implies that, in practice, a minor misspecification of the model for *G*(·|**Z**) can be corrected by doing a good job in modeling *F*(·|**Z**), and vice versa.

Writing * _{i}* (

$$\begin{array}{c}{\mathbf{U}}_{DR}(\mathit{\beta})=\sum _{i=1}^{n}\left[\frac{{\mathrm{\Delta}}_{i}{\stackrel{\sim}{\mathbf{D}}}_{i}(\beta ;{T}_{i},{\mathbf{Z}}_{i})}{\widehat{G}({X}_{i})}-{\int}_{0}^{\tau}{\widehat{\mathbf{Q}}}_{F}(t,{\mathbf{Z}}_{i})\frac{d{\widehat{M}}_{i}^{c}(t)}{\widehat{G}(t)}\right],\\ {\widehat{\mathbf{Q}}}_{F}(t,{\mathbf{Z}}_{i})=\frac{1}{\widehat{F}(t\mid {\mathbf{Z}}_{i})}{\int}_{t}^{\tau}{\stackrel{\sim}{\mathbf{D}}}_{i}(\mathit{\beta};u,{\mathbf{Z}}_{i})d\widehat{F}(u\mid {\mathbf{Z}}_{i}),\end{array}$$

(22)

and (*t*|**Z**) is an estimator of *F*(*t*|**Z**). This may seem like circular reasoning, because if we had methods to estimate *F*(*t*|**Z**), we would not need to develop this whole theory in the first place. The rationale here is that we can propose a simple working model for estimating *F*(*t*|**Z**). For example, we may use a parametric model or a Cox model even if they are not correct. An interesting fact is that we can even use model (5) for this purpose. That is, we try to obtain an estimator of *F*(*t*|**Z**) under model (5) (note that this estimator cannot be based on the censoring distribution assumption), plug it back to (22), and construct more efficient and robust estimating equations. This way, the new estimator will be consistent and asymptotically normal if either model (5) (or whatever models we may wish to posit as a working model for approximating *F*(*t*|**Z**)) or the censoring mechanism assumption is correct.

In this section we report some results from simulation studies conducted to assess the performance of the estimation procedures developed in Sections 2, 3, and 4 with the focus on estimation of the scalar *β*_{0}. For cases in which the censoring time *C* is independent of *T* and *Z*, we let *Z* ~ Bernoulli{0, 1} with probability 0.5 and *C* follow a uniform distribution [0, *V*], with *V* varying to yield desiring censoring percentages. For cases in which the censoring time depends on the covariate *Z*, *C* is generated from the Cox model (8) with *γ*_{0} = 1 and *λ*_{0}(*t*) being various constants for generating different censoring percentages. In both cases, care is taken to ensure the support of *C* is longer than the support of *T* for all *Z*. We considered several transformation functions *g,* but for the purpose of comparison we will only report in this section the results from *g*(*t*) = exp(*t*), which corresponds to the proportional mean residual life model (3) and *g*(*t*) = *t*, which corresponds to the additive mean residual life model (4). Data analysis results based on other choices of *g* can be found in Section 6. To generate *T* according to model (5), the true *β*_{0} is chosen to be 0 or 0.5 and the baseline function *m*_{0}(*t*) is taken from the Hall-Wellner family—in other words, *m*_{0}(*t*) = *g*^{−1}{(*D*_{1}*t* + *D*_{2})^{+}}, where *D*_{1} > − 1, *D*_{2} > 0, and *d*^{+} denotes *dI*(*d* ≥ 0) for any quantity *d*. In our simulation, we consider two scenarios. One is that *D*_{1} = − 0.5 and *D*_{2} = 0.5, which gives a uniform distribution [0, 1] when *Z* = 0. The other is *D*_{1} = −1/3 and *D*_{2} = 1, which corresponds to the survival function *F*(*t*) = (1 − *t*/3)^{2}, 0 ≤ *t* ≤ 3 when *Z* = 0 and mimics the exponential distribution resulting from the right skewness. When *Z* = 0.5, one can use the formula

$$F(t\mid Z)={\scriptstyle \frac{m(0\mid Z)}{m(t\mid Z)}}exp\left\{-{\int}_{0}^{t}{\scriptstyle \frac{du}{m(u\mid Z)}}\right\}$$

to compute the conditional survival function of *T*, where *m*(*t*|*Z*) follows model (5) and depends on the choice of *g*.

Because it is difficult to evaluate the effect of the integrator *H* analytically, we compare simulation results from different *H* with the focus on the variance of * _{a}* and

Table 1 shows the results from independent censoring. The censoring percentages are around 10% or 30% and the sample sizes are 100 or 200. It can be seen that the point estimates of *β*_{0} and their standard error estimates are very accurate. Besides, the 95% empirical coverage probabilities based on normal approximation are reasonable. Table 2 shows similar results for cases in which the censoring time depends on the covariate. For both cases, the baseline MRLF with *D*_{1} = −0.5 and *D*_{2} = 0.5 mentioned earlier is used. Q-Q plots for the estimates indicate that asymptotic normality is a good approximation (figures not shown). Therefore, we conclude that the limiting distributions for the proposed estimators work quite well with moderate sample sizes. We also considered other choices of the transformation function *g* and different simulation setups, and obtained similar results.

In Table 3 we compare the performance of IPCW and two AIPCW estimators using *H*_{1}(*t*). Moreover, we program the quasi-partial score (QPS) estimators proposed in Chen and Cheng (2005, 2006) for the proportional mean residual life model and the additive mean residual life model, respectively. To make our simulation results comparable with those of Chen and Cheng (2005, 2006), we mimic their simulation setup by letting the censoring be an independent uniform random variable, covariate a Bernoulli random variable as before, and the baseline *T* either uniform or right skewed as mentioned earlier. Because the validity of the IPCW and QPS estimators have already been established earlier and in Chen and Cheng (2005, 2006), respectively, we only report their relative efficiency with respect to the AIPCW estimator (AIPCW1, to be defined shortly). Here, by convention, the relative efficiency of an unbiased estimator _{1} with respect to another unbiased estimator _{2} for the true parameter *θ*_{0} is defined as var(_{2})/var(_{1}), which in our simulations is approximated by their corresponding sample estimates. In Table 3, the two AIPCW estimators are based on (22). That is, they are the efficient and DR estimators. The difference between them is that one (denoted by AIPCW1 in Table 3) uses the Kaplan-Meier estimator for calculating (*t*|*Z*), which is applicable when the covariate takes finitely many values. The other (denoted by AIPCW2 in Table 3) uses the estimate * _{a}*,

$$\begin{array}{l}\widehat{F}(t\mid Z)=\frac{\widehat{m}(0\mid Z)}{\widehat{m}(t\mid Z)}exp\left\{-{\int}_{0}^{t}\frac{du}{\widehat{m}(u\mid Z)}\right\}\\ =\frac{g\{{\widehat{m}}_{a0}(0)+{\widehat{\beta}}_{a}Z\}}{g\{{\widehat{m}}_{a0}(t)+{\widehat{\beta}}_{a}Z\}}exp\left\{-{\int}_{0}^{t}\frac{du}{g\{{\widehat{m}}_{a0}(u)+{\widehat{\beta}}_{a}Z\}}\right\}.\end{array}$$

Note that this would be the choice when the covariate is continuous or takes too many values compared with the sample size. It can be seen that both AIPCW estimators are quite accurate, with AIPCW1 slightly more efficient. Moreover, the AIPCW estimators are more efficient than the IPCW estimator as expected. The comparisons between our estimators and the QPS estimators indicate that our estimators have greater efficiency in all scenarios.

In this section, we apply the proposed methodology to two real-life datasets obtained from clinical trials.

The purpose of this clinical trial was to investigate the impact of the systematic combination of chemotherapy with regular radiotherapy on survival and recurrence patterns in patients with surgically resected lung cancer. In total, 172 patients were randomized between November 1979 and May 1985 with 86 in the radiotherapy (RT) group, 78 in the radio-therapy plus chemotherapy (RT+CAP) group, and the remaining eight patients became ineligible for statistical analysis. Both time to death and disease-free survival time are of interest in this study. For the purpose of illustration, we consider only the disease-free survival time as the response (*T*), which is subject to right censoring if the patient was alive and had never had a tumor recurrence by the end of the study. In addition to the treatment factor, which is denoted by *Z*_{1}, we also consider another important histology factor, squamous versus non-squamous/mixed cell type, and denote it by *Z*_{2}. We define *Z*_{1} = 1 if the patient was in the RT+CAP group and 0 otherwise, and *Z*_{2} = 1 if the patient had the squamous cell type and 0 otherwise. For a detailed description and some preliminary statistical analysis results, see Lad, Rubinstein, and Sadeghi (1988) and Piantadosi (2005).

To illustrate our method, we first apply the Cox regression model by using the censoring time as the response and the two covariates mentioned earlier. It appears that the covariate effects are not significant (both *p*-values > 0.50). However, we caution here that it is certainly possible to apply other approaches to find better solutions for modeling the censoring mechanism. Besides, it could be the case that the covariate effects (on censoring) are nonlinear or the distributions are nonnormal so that the regression coefficients are zero but independence fails, as pointed out by the editor. However, for purposes of illustration and for the time being, we treat censoring as independent for simplicity. We then apply the proposed transformed mean residual life model and estimate the regression parameter *β*_{0} (together with the baseline MRLF *g*{*m*_{0}(*t*)}) by using the estimating equations (6) and (7) with the integrator *H*_{1} and three different choices of the transformation function *g*. It is seen (Table 4) that both covariates have a significant impact on the mean residual life. Namely, with adjuvant chemotherapy, the expected residual disease-free time is significantly longer. Also, patients in the nonsquamous/mixed cell-type group tend to have longer expected residual disease-free time.

Notice that the magnitudes of the parameter estimates are quite different for various choices of the link function *g*. This is because the parameters have different interpretations for different *g*. For instance, if *g*(*t*) = *t* (i.e., model (5) reduces to the additive model (4)), then the parameter estimate means that the adjuvant chemotherapy increases the expected residual disease-free time by approximately 383 days when compared with radiotherapy alone. On the other hand, when *g*(*t*) = exp(*t*), which yields model (3), the parameter estimate should be interpreted as the ratio of the mean residual function between the two groups. In other words, the adjuvant chemotherapy doubles (*e*^{0.726} ≈ 2) the mean residual disease-free time when compared with radiotherapy alone. Also using the proportional mean residual life model and their estimation approach, Chen et al. (2005) obtained similar results.

An interesting phenomenon is that if we use the Cox proportional hazards model for this problem we will get a different interpretation on the effect of cell type (although for the effect of treatment, the interpretations are in the same direction). The Cox model indicates that both treatment and cell type are significant covariates (*p*-values are 0.001 and 0.003, respectively), but the coefficient for cell type is −0.570, which means patients in the nonsquamous/mixed cell-type group have a higher risk. A closer examination of the data reveals the reason. The largest uncensored observation (recurrence time) in the squamous group is only 1,388 days, whereas the largest three uncensored observations in the nonsquamous/mixed group are 1,577, 2,375, and 2,566 days. In other words, the mean residual time for the nonsquamous/mixed group is greatly inflated by these extremely large values. One can also see this from the fact that the Cox proportional hazards model and most other survival regression models based on the hazard function essentially measure the instantaneous risk for failure, which might be explained as a conditional probability, whereas the mean residual lifetime measures some conditional average life span and thus can be heavily affected by extreme values. As is well known, the Cox model parameter estimate is rank invariant, but that of the mean residual lifetime model is not.

Theoretically, it has been pointed out (Shaked and Shantikhumar, 1993) that the ordering of the MRLF is weaker than the ordering of the hazard function. Namely, if two MRLFs *m _{a}*(

The interest of a recent study conducted at the Memorial Sloan-Kettering Cancer Center is to examine the effect of the radiation therapy dose level and patient age on the PSA relapse time (see the Introduction for the definition). Let *Z*_{1} denote the age of the patient before the initial treatment, which ranges from 42 to 87 years continuously, and *Z*_{2} denote dose, which takes two values: 1 = low dose (≤70.2 Gy) and 2 = high dose ((75.6 Gy, 81 Gy]) (these cutoffs are suggested by experienced radiotherapists). In total, 191 patients underwent radiotherapy. These patients were all classified as high risk based on their baseline PSA value, Gleason score, and T-stage. They also all had adjuvant hormonal therapy. In this analysis, we first apply a Cox model on the censoring time using age and dose level as covariates. It turns out that both are significant and thus we have the situation discussed in Section 3 (i.e., the censoring mechanism depends on covariates and needs to be modeled). Again, we emphasize here that a more comprehensive study may be conducted to fine-tune the modeling approaches for the censoring distribution. However, for the time being, we shall use the Cox model for such a purpose. Implementing (10) and (11) gives us coefficient estimates of age and dose level of 14.60 and 127.47, respectively. The standard errors for these two point estimates, based on Theorem 2, are 8.61 and 68.12. These results suggest that higher doses significantly elongate the residual PSA relapse time, whereas age has only a marginal effect. An interesting issue is that older patients tend to have longer relapse times, which is somewhat counterintuitive, but agrees with a previous study using the Cox model for assessing hazards (Zelefsky et al. 2008). Note here that, as discussed in the Introduction, we used the Box-Cox transformation for the link function *g* with *ρ* = 0.5 because neither the additive nor the multiplicative modeling assumption sounds natural for this clinical problem. However, one can certainly enjoy other choices of *g*, and thus the flexibility of our approach, based on prior data information or physicians’ assertions if they are available.

In this article we proposed a class of transformation mean residual life models. This very general class of models encompasses some previously established models. For example, the proportional mean residual life model and the additive mean residual life model are special cases of our model. IPCW-based generalized estimating equation approaches are used to estimate the regression parameters of interest. We also considered efficiency and double robustness properties of our estimators and suggested improved (AIPCW) estimators. Both asymptotic theory and simulation studies show that the IPCW and the AIPCW estimators perform well.

Notice that, unlike some previously established methods, the proposed estimating equations here have intuitive meanings. They essentially yield the weighted least-square-type estimators, with the weight determined by the censoring probability. Therefore, although our approaches are proposed for regression purposes, they can be very easily modified for one-sample problems. Our approaches provide a nice platform for obtaining efficient estimators (without censoring, our approaches naturally reduce to that of Bickel et al. (1993)), and our efficiency and double robustness arguments also seem to be the first for the MRLF estimation problems under right censoring.

As we saw in both theory and applications, understanding and modeling the censoring mechanism is important for our approaches. Even though the double robustness property provides us a certain amount of luxury of misspecifying the model of censorship, we should make every effort to find a good approximation. In our numerical studies, when the censoring time may depend on **Z**, the Cox proportional hazards model was used as the working model for the censoring time. Of course, we can also choose some other semiparametric regression models as the working model for censoring. For example, we may use the proportional mean residual life model or the additive mean residual life model and obtain the estimators of the censoring model parameters using the approach of Chen and Cheng (2005) or Chen and Cheng (2006). Then the estimator of the censoring survival distribution *G*(*t*; **Z**) can be obtained using the inversion formula as before. Thus, the unknown parameter in model (5) can be estimated by using the procedures in Sections 3 and 4.

For the examples, various models with different choices of the link function *g* were considered. It would be helpful to provide a method for selecting or comparing different link functions. For a given link function, following Lin et al. (2000), we can use a residual-based procedure for checking the adequacy of the model. If, on the other hand, we restrict the function *g* to the Box-Cox transformation family but allow *ρ* to be estimated from the observed data, then it would provide a way in model selection between the additive mean residual life model and other alternatives. Let *g _{ρ}*(

$$\begin{array}{l}{\mathbf{U}}_{a}^{\ast}(\mathit{\beta},\rho )=\sum _{i=1}^{n}{\int}_{0}^{\tau}\frac{{\mathrm{\Delta}}_{i}I({X}_{i}>t){\mathbf{Z}}_{i}^{\ast}(t)}{\widehat{G}({X}_{i})}[({X}_{i}-t)\\ -{g}_{\rho}\{{\widehat{m}}_{a0}^{\ast}(t;\mathit{\beta},\rho )+{\beta}^{\prime}{\mathbf{Z}}_{i}\}]dH(t)=\mathbf{0}\end{array}$$

under independent censoring and

$$\begin{array}{l}{\mathbf{U}}_{b}^{\ast}(\mathit{\beta},\rho )=\sum _{i=1}^{n}{\int}_{0}^{\tau}\frac{{\mathrm{\Delta}}_{i}I({X}_{i}>t){\mathbf{Z}}_{i}^{\ast}(t)}{{G}_{i}({X}_{i};\widehat{\mathit{\gamma}},{\widehat{\mathrm{\Lambda}}}_{0})}[({X}_{i}-t)-{g}_{\rho}\{{\widehat{m}}_{b0}^{\ast}(t;\mathit{\beta},\rho )\\ +{\mathit{\beta}}^{\prime}{\mathbf{Z}}_{i}\}]dH(t)=\mathbf{0}\end{array}$$

under dependent censoring, where
${\mathbf{Z}}_{i}^{\ast}(t)=({{\mathbf{Z}}^{\prime}}_{i},{W}_{i}(t){)}^{\prime}$, *W _{i}*(

One important issue with the mean residual life model is the embedded model constraint. That is, *m*(*t*) + *t* should be nondecreasing in *t*. Chen and Cheng (2005, 2006) discussed this constraint for the proportional mean residual life model (3) and the additive mean residual life model (4). For the class of transformed mean residual life model (5), however, *m*(*t*|**Z**) + *t* = *g*{*m*_{0}(*t*) + ** β**′

Another interesting issue is the effect of the integrator *H*, which essentially plays the role of a weight function. Our simulation results seem to suggest that estimating the baseline function *m*_{0}(*t*) at all uncensored observations does not necessarily improve either the accuracy of the point estimate or its efficiency. Here we refer interested readers to Lin, Wei, and Ying (2001), where they also pointed out the difficulty of finding an optimal *H* analytically and thus compared different *H* by simulation in dealing with point processes.

The first author’s research was supported in part by the National Natural Science Foundation of China Grants (nos. 10571169 and 10731010) and the National Basic Research Programme (no. 2007CB814902). The second author’s research was supported in part by the National Cancer Institute Core Grant (CA008748). The authors thank the editor, Professor Stephen L. Portnoy, an associate editor, and two referees for their insightful comments and suggestions that greatly improved the article.

The proofs of asymptotic properties and double robustness are provided in the supplemental technical report. Here, we only outline the regularity conditions.

- (C1)
*H*(*t*) converges almost surely to a nonrandom and bounded function (*t*) uniformly in*t*[0,*τ*]. - (C2)
**Z**is bounded and*G*is continuous. - (C3)
*m*_{0}(*t*) is continuously differentiable on [0,*τ*]. - (C4) $\mathbf{A}:=E[{\int}_{0}^{\tau}I({T}_{i}>t){\{{\mathbf{Z}}_{i}-{\overline{\mathbf{z}}}_{a}(t)\}}^{\otimes 2}\stackrel{.}{g}\{{m}_{0}(t)+{\beta}_{0}^{\prime}{\mathbf{Z}}_{i}\}d\stackrel{\sim}{H}(t)]$is nonsingular where
(_{a}*t*) is the limit of(_{a}*t*;*β*_{o}). - (C5) $\mathbf{B}:=E[{\int}_{0}^{\tau}I({T}_{i}>t){\{{\mathbf{Z}}_{i}-{\overline{\mathbf{z}}}_{b}(t)\}}^{\otimes 2}\stackrel{.}{g}\{{m}_{0}(t)+{\beta}_{0}^{\prime}{\mathbf{Z}}_{i}\}d\stackrel{\sim}{H}(t)]$is nonsingular where
(_{b}*t*) is the limit of(_{b}*t*;*β*_{o}).

Liuquan Sun, Professor, Institute of Applied Mathematics, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100080, P.R. China (E-mail: nc.ca.tma@qls)

Zhigang Zhang, Assistant Attending Biostatistician, Department of Epidemiology and Biostatistics, Memorial Sloan-Kettering Cancer Center, New York, NY 10065 (E-mail: gro.ccksm@zgnahz)

- Arnold BC, Zahedi H. On Multivariate Mean Remaining Life Functions. Journal of Multivariate Analysis. 1988;25:1–9.
- Balkema AA, de Haan L. Residual Life Time at Great Age. Annals of Probability. 1974;2:792–804.
- Bang H, Tsiatis AA. Estimating Medical Costs with Censored Data. Biometrika. 2000;87:329–343.
- Bang H, Tsiatis AA. Median Regression with Censored Cost Data. Biometrics. 2002;58:643–649. [PubMed]
- Bassan B, Kochar S, Spizzichino F. Some Bivariate Notions of IFR and DMRL and Related Properties. Journal of Applied Probability. 2002;39:533–544.
- Bickel PJ, Klaassen CAJ, Ritov Y, Wellner JA. Efficient and Adaptive Inference in Semiparametric Models. Baltimore, MD: Johns Hopkins University Press; 1993.
- Chen YQ. Additive Regression of Expectancy. Journal of the American Statistical Association. 2007;102:153–166.
- Chen YQ, Cheng SC. Semiparametric Regression Analysis of Mean Residual Life with Censored Survival Data. Biometrika. 2005;92:19–29.
- Chen YQ, Cheng SC. Linear Life Expectancy Regression with Censored Data. Biometrika. 2006;93:303–313.
- Chen YQ, Jewell NP, Lei X, Cheng SC. Semiparametric Estimation of Proportional Mean Residual Life Model in Presence of Censoring. Biometrics. 2005;61:170–178. [PubMed]
- Cox DR. Renewal Theory. London: Methuen; 1961.
- Dauxois J-Y. Estimating the Deviation from Exponentiality under Random Censorship. Communications in Statistics Theory and Methods. 2003;32:2117–2137.
- El-Bassiouny AH, Alwasel IA. A Goodness of Fit Approach to Testing Mean Residual Times. Applied Mathematics and Computation. 2003;143:301–307.
- Fleming TR, Harrington DP. Counting Processes and Survival Analysis. New York: John Wiley; 1991.
- Grambsch P, Therneau T. Proportional Hazards Tests and Diagnostics Based on Weighted Residuals. Biometrika. 1994;81:515–526.
- Hall WJ, Wellner JA. Mean Residual Life. In: Csörgő M, Dawson DA, Rao JNK, Saleh AKME, editors. Proceedings of the International Symposium on Statistics and Related Topics. Amsterdam: North Holland; 1984. pp. 169–184.
- Hollander M, Proschan F. Tests for the Mean Residual Life. Biometrika. 1975;62:585–593.
- Kotz S, Shanbhag DN. Some New Approaches to Probability Distribution. Advances in Applied Probability. 1980;12:903–921.
- Lad T, Rubinstein L, Sadeghi A. The Benefit of Adjuvant Treatment of Resected Locally Advance Non-Small-Cell Lung Cancer. Journal of Clinical Oncology. 1988;6:9–17. [PubMed]
- Liang KY, Zeger SL. Longitudinal Data Analysis Using Generalized Linear Models. Biometrika. 1986;73:13–22.
- Lin DY, Wei LJ, Yang I, Ying Z. Semiparametric Regression for the Mean and Rate Functions of Recurrent Events. Journal of the Royal Statistical Society, Ser B. 2000;62:711–730.
- Lin DY, Wei LJ, Ying Z. Semiparametric Transformation Models for Point Processes. Journal of the American Statistical Association. 2001;96:620–628.
- Maguluri G, Zhang CH. Estimation in the Mean Residual Life Regression Model. Journal of the Royal Statistical Society, Ser B. 1994;56:477–489.
- Newey WK. Semiparametric Efficiency Bounds. Journal of Applied Econometrics. 1990;5:99–135.
- Oakes D, Dasu T. A Note on Residual Life. Biometrika. 1990;77:409–410.
- Piantadosi S. Clinical Trials: A Methodologic Perspective. 2. New York: Wiley; 2005.
- Robins JM, Ritov Y. Towards a Curse of Dimensionality Appropriate (CODA) Asymptotic Theory for Semiparametric Models. Statistics in Medicine. 1997;16:285–319. [PubMed]
- Robins JM, Rotnitzky A. Recovery of Information and Adjustment for Dependent Censoring Using Surrogate Markers. In: Jewell N, Dietz K, Farewell V, editors. AIDS Epidemiology: Methodologic Issues. Boston: Birkhauser; 1992. pp. 297–331.
- Robins JM, Rotnitzky A. Comment on the Bickel and Kwon Article, ‘Inference for Semiparametric Models: Some Questions and an Answer’ Statistica Sinica. 2001;11:920–936.
- Robins JM, Rotnitzky A, Zhao LP. Estimation of Regression Coefficients When Some Regressors Are Not Always Observed. Journal of the American Statistical Association. 1994;89:846–866.
- Rojo J, Ghebremichael M. Estimation of Two Ordered Bivariate Mean Residual Life Functions. Journal of Multivariate Analysis. 2006;97:431–454.
- Rotnitzky A, Robins JM. Inverse Probability Weighting in Survival Analysis. In: Armitage P, Colton T, editors. Encyclopedia of Biostatistics. 2. Chichester, UK: Wiley; 2005. pp. 2619–2625.
- Rubin D, van der Laan MJ. Doubly Robust Censoring Unbiased Transformations. 2006. available at http://works.bepress.com/mark_van_der_laan/57. [PubMed]
- Scharfstein DO, Rotnitzky A, Robins JM. Adjusting for Nonignorable Drop-Out Using Semiparametric Nonresponse Models,” (with discussions and rejoinder) Journal of the American Statistical Association. 1999;94:1096–1146.
- Shaked M, Shantikhumar J. Stochastic Orders and Their Applications. Boston: Academic Press; 1993.
- Swartz GB. The Mean Residual Life Function. IEEE Transactions on Reliability. 1973;R-22:108–109.
- Tsiatis AA. Semiparametric Theory and Missing Data. New York: Springer-Verlag; 2006.
- van der Laan MJ, Robins JM. Unified Methods for Censored Longitudinal Data and Causality. New York: Springer Verlag; 2003.
- Zelefsky MJ, Yamada Y, Fuks Z, Zhang Z, Hunt M, Cahlon O, Park J, Shippy A. Long-Term Results of Conformal Radiotherapy for Prostate Cancer: Impact of Dose Escalation on Biochemical Tumor Control and Distant Metastases-Free Survival Outcomes. International Journal of Radiation Oncology Biology Physics. 2008;71:1028–1033. [PubMed]
- Zhao H, Tsiatis AA. Efficient Estimation of the Distribution of Quality-Adjusted Survival Time. Biometrics. 1999;55:1101–1107. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |