Efficient estimation of accelerated lifetime models under length-biased sampling

TL;DR

This paper investigates efficient methods for estimating accelerated lifetime models in length-biased sampling scenarios, demonstrating that naive estimators can be fully efficient when covariate distribution is unspecified.

Contribution

It shows that in length-biased sampling, naive estimators based on conditional likelihood are fully efficient if the covariate distribution is unspecified.

Findings

Naive estimators are fully efficient under certain conditions.

The invariance of the AFT model under length-biased sampling.

Efficiency depends on covariate distribution assumptions.

Abstract

In prevalent cohort studies where subjects are recruited at a cross-section, the time to an event may be subject to length-biased sampling, with the observed data being either the forward recurrence time, or the backward recurrence time, or their sum. In the regression setting, it has been shown that the accelerated failure time model for the underlying event time is invariant under these observed data set-ups and can be fitted using standard methodology for accelerated failure time model estimation, ignoring the length-bias. However, the efficiency of these estimators is unclear, owing to the fact that the observed covariate distribution, which is also length-biased, may contain information about the regression parameter in the accelerated life model. We demonstrate that if the true covariate distribution is completely unspecified, then the naive estimator based on the conditional likelihood given the covariates is fully efficient.

Tables

Table 1. Table 1: Estimates for the Backward Recurrence Time Data

Parameters	Sample	Profile Likelihood Approach			Vanishing Mean		Known Covariance
	Size	Bias	SE	CP (%)	Bias	SE	Bias	SE
$\theta =1$	100	$-0.033$	0.258	94.4	$-0.110$	0.230	$-0.051$	0.201
	200	0.006	0.188	95.6	$-0.052$	0.161	$-0.001$	0.137
	400	0.010	0.133	95.1	0.002	0.111	0.001	0.099
$\theta =0.5$	100	0.028	0.202	95.8	$-0.165$	0.190	$-0.020$	0.142
	200	0.016	0.188	95.1	$-0.033$	0.134	0.007	0.095
	400	0.010	0.117	95.5	0.015	0.092	0.004	0.067
$\theta =2$	100	0.010	0.299	93.9	0.027	0.553	0.027	0.366
	200	0.003	0.203	95.7	0.055	0.284	0.004	0.216
	400	$-0.005$	0.169	94.5	0.012	0.198	$-0.003$	0.105
$x=-0.9$	100	$-0.029$	0.254	95.2	$-0.055$	0.250	0.031	0.222
$\theta =1$	200	0.009	0.181	94.1	0.181	0.171	0.115	0.138
	400	0.002	0.124	94.6	0.155	0.114	0.213	0.094
$x=-0.8$	100	0.023	0.269	93.3	0.335	0.270	0.195	0.211
$\theta =1$	200	$-0.006$	0.169	94.9	0.490	0.184	0.234	0.142
	400	0.003	0.108	94.4	0.304	0.127	0.278	0.099

Table 2. Table 2: Estimates for time ratios and the corresponding confidence intervals for nulliparous women

		Semiparametric AFT	OLS
Covariate	No	Time Ratio	Time Ratio
Tobacco Consumption at recruitment
Non-Smokers	159	1	1
Smokers	92	1.20(0.75,1.78)	1.04(0.70,1.53)
Age at recruitment
0-17	3	7.50(1.50,38.0)	7.32(1.29,41.4)
18-24	50	2.00(1.20,3.41)	2.08(1.24,3.49)
25-29	93	1	1
30-34	62	1.00(0.61,1.74)	1.01(0.63,1.64)
35-39	41	1.10(0.61,2.02)	0.93(0.54,1.62)
40-44	2	0.13(0.01,1.17)	0.13(0.02,1.10)
Frequency of Sexual Intercourse
1 per month	0
1-3 per month	44	2.20(1.20,3.89)	2.18(1.27,3.71)
1-2 per week	109	1.20(0.78,1.92)	1.23(0.81,1.86)
$\ge$3 per week	98	1	1
Menstrual Cycle Length
27 days	53	1	1
27-29 days	110	0.90(0.52,1.55)	0.86(0.52,1.41)
$\ge$30 days	88	1.10(0.62,1.81)	0.99(0.59,1.63)

Equations

$$g_{\tilde{T}}(t)=S_T(t){\mu }_T^{-1},t\in (0,\infty ),$$

$$f_{T_f|T_{\mathrm{LB}}}(s|t)=I(0<s<t)t^{-1},$$

$$T=e^{{\theta }^{\prime }Z}U,$$

$$f_{T_{LB},\tilde{Z}}(t,z)=t{e}^{-{\theta }^{\prime }z}g(e^{-{\theta }^{\prime }z}t){\mu }_g^{-1}h(z)/\int e^{{\theta }^{\prime }z}h(z)dz,$$

$$T_{\mathrm{LB}}=e^{{\theta }^{\prime }{Z}_{\mathrm{LB}}}{U}_{\mathrm{LB}},$$

$$h_{Z_{\mathrm{LB}}}(z)=e^{{\theta }^{\prime }z}h(z)/\int e^{{\theta }^{\prime }z}h(z)dz$$

$$f_{\tilde{T},\tilde{Z}}(t,z)=\frac{e^{-{\theta }^{\prime }z}S(e^{-{\theta }^{\prime }z}t)}{{\mu }_g}\times \frac{e^{{\theta }^{\prime }z}h(z)}{\int e^{{\theta }^{\prime }u}h(u)du}.$$

$$\tilde{T}=e^{{\theta }^{\prime }\tilde{Z}}\tilde{U},$$

$${\mathcal{P}}^{\ast }=\left\{P_{\theta ,g,h}^{\ast }:\theta \in \Theta ,g\in \mathcal{G},h\in \mathcal{H}\right\},$$

$$\frac{d{P}_{\theta ,g,h}^{\ast }}{d\nu }(t,z)=e^{(\theta -{\theta }_0{)}^{\prime }z}g\{e^{(\theta -{\theta }_0{)}^{\prime }z}t\}h(z),$$

$$\mathcal{P}=\left\{P_{\theta ,S_g,h}:\theta \in \Theta ,g\in \mathcal{G},h\in {\mathcal{H}}^{\prime }\right\},$$

$${\mathcal{H}}^{\prime }=\left\{h:h\in \mathcal{H},\int e^{{\theta }^{\prime }z}h(z)d(z)<\infty ,\int z^{\otimes 2}{e}^{{\theta }^{\prime }z}h(z)dz<\infty ,\theta \in \Theta \right\}.$$

$$\frac{d{P}_{\theta ,S_g,h}}{d\nu }(t,z)=\frac{e^{(\theta -{\theta }_0{)}^{\prime }z}{S}_g\{e^{(\theta -{\theta }_0{)}^{\prime }z}t\}}{\int S_g(v)dv}\times \frac{e^{{\theta }^{\prime }z}h(z)}{\int e^{{\theta }^{\prime }z}h(z)dz}.$$

$$g_{U(\theta )}(u)=\frac{e^{(\theta -{\theta }_0{)}^{\prime }z}S\{e^{(\theta -{\theta }_0{)}^{\prime }z}u\}}{\int S(v)dv},$$

$${\lambda }_{U(\theta )}(u)=\frac{S\{e^{(\theta -{\theta }_0{)}^{\prime }z}u\}}{{\int }_u^{\infty }S\{e^{(\theta -{\theta }_0{)}^{\prime }z}w\}dw},u\in (0,\infty ).$$

$$M(t)=I\{U(\theta )\le t\}-{\int }_0^{t}I\{U(\theta )>s\}{\lambda }_{U(\theta )}(s)ds$$

$$Ra(t)=a(t)-\frac{{\int }_t^{\infty }a(u)S(u)du}{{\int }_t^{\infty }S(u)du},\text{ for }a\in L_2^0(S).$$

$${\dot{l}}_{{\theta }_0}=\tilde{Z}{\int }_0^{U^c({\theta }_0)}R\varphi (s)dM(s)+(\tilde{Z}-E\tilde{Z}),$$

$${\dot{l}}_{S}b={\int }_0^{U^c({\theta }_0)}Rb(s)dM(s),$$

$${\tilde{l}}_{\theta ,S}={\int }_0^{U^c({\theta }_0)}[\tilde{Z}-E\{\tilde{Z}|U^c({\theta }_0)\ge s\}]R\varphi (s)dM(s).$$

$$l(\theta )={\left\{g_{U(\theta )}(U_i)\right\}}^{{\delta }_i}{\left\{{\int }_{U_i^c}^{\infty }{g}_{U(\theta )}(u)du\right\}}^{1-{\delta }_i}{h}_{Z,\theta }({\tilde{Z}}_i).$$

$${\dot{l}}_{{\theta }_0}={\tilde{Z}}_i\left[{\delta }_i\varphi (U_i)+(1-{\delta }_i)E\left\{\varphi (Ut)\mid U(\theta )>U_i^c\right\}\right]+({\tilde{Z}}_i-E\tilde{Z}).$$

$$\eta \mapsto {\lambda }_{\eta }(t)=\frac{S_{\eta }(t)}{{\int }_t^{\infty }{S}_{\eta }(v)dv}.$$

$$E\left\{\left({\dot{l}}_{{\theta }_0}-{\dot{l}}_S{b}^{\ast }\right){\dot{l}}_{S}b\right\}$$

$$=EE\left\{({\dot{l}}_{{\theta }_0}-{\dot{l}}_S{b}^{\ast }){\dot{l}}_{S}b\mid \tilde{Z},U^c({\theta }_0)\right\}$$

$$=EE\left\{{\int }_0^{U^c({\theta }_0)}(\tilde{Z}R\varphi (s)-R{b}^{\ast }(s))Rb(s)I\{U({\theta }_0)\ge s\}{\lambda }_{U({\theta }_0)}(s)ds\mid \tilde{Z},U^c({\theta }_0)\right\}$$

$$=E\left\{{\int }_0^{U^c({\theta }_0)}(\tilde{Z}R\varphi (s)-R{b}^{\ast }(s))Rb(s)d{F}_{U({\theta }_0)}(s)\right\}$$

$$=\int \left\{E(\tilde{Z}I\{U^c({\theta }_0)\ge s\})R\varphi (s)-EI\{U^c({\theta }_0)\ge s\}R{b}^{\ast }(s)\right\}Rb(s)d{F}_{U({\theta }_0)}(s).$$

$$E{Y}_1{Y}_2=E\int f_1{f}_{2}d⟨M,M⟩=E\int f_1(s)f_2(s)I\{U({\theta }_0)\ge s\}d\Lambda (s).$$

$$R{b}^{\ast }(s)=\frac{E\{\tilde{Z}I(U^c({\theta }_0)\ge s)\}}{EI(U^c({\theta }_0)\ge s)}R\varphi (s)=E(\tilde{Z}|U^c({\theta }_0)\ge s)R\varphi (s).$$

$${\Pi }_0({\dot{l}}_{{\theta }_0}|{\dot{\mathcal{P}}}_S)={\int }_0^{U^c({\theta }_0)}E(\tilde{Z}|U^c({\theta }_0)\ge s)R\varphi (s)dM(s).$$

$${\dot{l}}_{h}k=k-\frac{\int e^{{\theta }_0^{\prime }z}k(z)h(z)dz}{\int e^{{\theta }_0^{\prime }z}h(z)dz}\equiv m(z),$$

Full text

accelerated failure time model,

backward recurrence time,

forward recurrence time,

length-biased time,

and

t1This work was done prior to the author joining the FDA and does not represent the official position of the FDA

1. 1 Introduction

2. 2 Efficient Scores and Estimation

3. 2.1 Inference for Forward and Backward Recurrence Times

4. 2.2 Inference for Length-Biased Data

5. 3 Simulation Studies

6. 4 Data Analysis

7. 5 Discussion

1 Introduction

Many study designs involve cross-sectional sampling, which may lead to length-biased sampling of a time to event $T$. In a prospective cohort study where the initiation time for the event is unknown and subjects are followed prospectively, a right-censored forward recurrence time $T_{f}$ is observed. This occurs, for example, in HIV seroprevalence studies [2], where time to AIDS following infection with HIV is of interest but the infection time is unknown. If the initiation time is known, but there is no follow-up, the backward recurrence time $T_{b}$ is observed. This current duration design [9] has been employed in pregnancy surveys, where current trier couples provide the length of an ongoing attempt at pregnancy, and in mover-stayer models [16]. If there is both an initiation time and a follow-up, a biased event time $T_{\mathrm{LB}}=T_{f}+T_{b}$ may be observed. When the sampling time is known, such data are commonly analyzed using methods for left-truncated data, where one conditions on the lack of an event prior to the sampling time, that is, $T_{\mathrm{LB}}>T_{b}$.

In the aforementioned data set-ups, only subjects who have experienced an initiating event prior to sampling can potentially be sampled, and the sample is biased towards larger values of $T$. If one assumes that the rate of the initiating event is stationary over time, e.g., is a homogeneous Poisson process, then the sampling time falls uniformly in the interval between the initiation and the event times [6, 15]. Letting $F_{T}$ denote the distribution of $T$, the length-biased version $T_{\mathrm{LB}}$ has distribution $F_{\mathrm{LB}}(t)=\int_{0}^{t}udF_{T}(u)\mu_{T}^{-1}$, $\quad t\geq 0$, where $\mu_{T}=\int_{0}^{\infty}udF_{T}(u)$. Under the uniform sampling time assumption [6, 14, 9], $T_{f}=T_{\mathrm{LB}}V$, where $V$ is uniform(0,1) and independent of $T_{\mathrm{LB}}$. Thus, both $T_{f}$ and $T_{b}$ will have the same density function. Since the density is the same, we use $\tilde{T}$ to denote both $T_{f}$ and $T_{b}$. The density of $\tilde{T}$ is

[TABLE]

where $S_{T}=1-F_{T}$ is the survival function of $T$. This well-known result, given in expression (2) of [8], can be derived from the uniformity of $V$, which yields that the conditional density of $T_{f}$ given $T_{\mathrm{LB}}$ has the form

[TABLE]

yielding the joint density $f_{T_{f},T_{\mathrm{LB}}}(s,t)=I(0<s<t)f_{T}(t)\mu_{T}^{-1}$, from which (1.1) follows.

In the presence of a $p\times 1$ covariate vector $Z$ with density $h$, one may formulate the effect of $Z$ on the underlying event time $T$ via the accelerated lifetime model

[TABLE]

where $\theta$ is a $p\times 1$ regression parameter and $U$ is a non-negative random variable with density $g$, survival function $S$ and hazard function $\lambda(u)=g(u)/S(u)$. In the semiparametric version of (1.3), with the distribution of $U$ completely unspecified, efficient estimation of $\theta$ without length-bias is achievable with and without right censoring [17].

The observed covariate is also subject to length-biased sampling. Using arguments in [4] and [13], we first obtain that a consequence of (1.3) is that the joint density of $T_{LB}$ and observed covariate $\tilde{Z}$ is proportional to $te^{-\theta^{\prime}z}g(e^{-\theta^{\prime}z}t)h(z)$ due to the length bias, yielding the joint density

[TABLE]

where the integral is over the range of $Z$ and $\mu_{g}$ is the mean associated with density $g$. This means that $T_{LB}$ satisfies an accelerated life model, given by

[TABLE]

where $Z_{\mathrm{LB}}$ is independent of $U_{\mathrm{LB}}$ with density of the form

[TABLE]

and $U_{\mathrm{LB}}$ has the possibly non-monotone density $g_{U_{\mathrm{LB}}}(u)=ug(u)\mu_{g}^{-1}$, $u\in(0,\infty)$.

The relation $\tilde{T}=T_{LB}V$ still holds for a uniform$(0,1)$ $V$ independent of $T_{LB}$, so we can multiply (1.4) by (1.2), integrate over $t$, and replace $s$ with $t$, to obtain the joint density of $\tilde{T}$ and observed covariate $\tilde{Z}$:

[TABLE]

Thus, if $T$ follows model (1.3), then the distribution of $\tilde{T}$ also follows an accelerated lifetime model

[TABLE]

where $\tilde{Z}$ has a density of the form $h_{\tilde{Z},\theta}(z)=e^{\theta^{\prime}z}h(z)/\int e^{\theta^{\prime}u}h(u)du$ and $\tilde{U}$ has monotone density $g_{\tilde{U}}(u)=S(u)/\int_{0}^{\infty}S(v)dv$. Thus the accelerated lifetime structure is maintained in models for the forward and backward recurrence times, as discussed in [10], as well as for the length-biased time.

Since the conditional distributions of $\tilde{T}$ and $T_{\mathrm{LB}}$ satisfy accelerated lifetime models, existing estimation procedures may naively be applied to obtain semiparametric estimators for $\theta$. However, as the marginal distribution of the observed covariates depends on the parameter $\theta$, it is unclear whether estimators derived from the conditional distributions of $\tilde{T}$ and $T_{\mathrm{LB}}$ will be fully efficient. Estimation using observed covariates has been considered in [5], and [3]. Under restrictive assumptions on $h$, for example, moment restrictions, improved estimation is possible. However, a comprehensive study of such issues with completely unspecified covariate distribution and right censoring has not been undertaken for general length-biased sampling.

The main contribution of this paper is to show that a naive efficient estimator which ignores the dependence on $\theta$ in the marginal covariate distribution is still efficient for estimation of the regression parameter in length-biased and recurrence time data. Hence, the standard techniques that are used for estimation in the accelerated failure time models can also be applied in these cases without loss of information. We provide the theoretical derivation of the efficient score in Section 2, simulation results and data analysis in Section 3 and 4 and conclude with some discussions in Section 5.

2 Efficient Scores and Estimation

We start by defining some of the important assumptions and notation used in the paper, along with a description of some of the tangent spaces used in deriving the efficient score. The concepts and notation closely follow that given in Chapters 3 and 18 of [12].

Let $\mathcal{G}$ be the class of all density functions on $\Re^{+}$ and $\mathcal{H}$ be the class of all density functions on $\Re^{p}$. The semiparametric model for the core accelerated lifetime model is given by

[TABLE]

where the distribution $P_{\theta,g,h}^{*}$ has a density with respect to an absolutely continuous measure $\nu$,

[TABLE]

where $\theta_{0}$ is the true parameter value. For the accelerated failure time model of the recurrence times, the semiparametric model is

[TABLE]

where $S_{g}(u)=\int_{u}^{\infty}g(v)dv$ for $g\in\mathcal{G}$, and

[TABLE]

We assume that $\Theta$ is a compact subset of $\mathcal{R}^{p}$. Further, we assume $P_{\theta,S_{g},h}$ has density

[TABLE]

Define $\mathcal{S}=\left\{S_{g}:g\in\mathcal{G}\right\}$. Let the true distribution be $P_{0}=P_{\theta_{0},S_{0},h_{0}}$ with $S_{0}=S_{g_{0}}$. Define the separate submodels for each parameter, holding the other parameters fixed, as $\mathcal{P}_{\theta}=\{P_{\theta,S_{0},h_{0}}:\theta\in\Theta\}$, $\mathcal{P}_{S}=\{P_{\theta_{0},S,h_{0}}:S\in\mathcal{S}\}$ and $\mathcal{P}_{h}=\{P_{\theta_{0},S_{0},h}:h\in\mathcal{H}^{\prime}\}$.

Let $\dot{\mathcal{P}}_{\theta}$, $\dot{\mathcal{P}}_{S}$ and $\dot{\mathcal{P}}_{h}$ be the tangent spaces for $\mathcal{P}_{\theta},\mathcal{P}_{S}\textrm{ and }\mathcal{P}_{h}$ at $P_{0}=P_{\theta_{0},S_{0},h_{0}}$. By definition of tangent spaces in Chapter 18 of [12], these are all closed subsets of $L_{2}^{0}(P_{0})$, where $L_{2}^{0}(P_{0})$ denotes square-integrable functions integrating to zero with respect to $P_{0}$. For a density $\nu$, let $L_{2}(\nu)$ denote the space of square-integrable functions with respect to the measure $\int\nu$. For a survival function $S$, we similarly let $L_{2}(S)$ denote the space of square-integrable functions with respect to the measure $\int S$, even though $S$ may not integrate to 1. Tangent spaces for a given model represent the set of all likelihood score functions for one-dimensional submodels of the given model. The three tangent spaces just defined represent the likelihood based scores used to estimate the parameter given in the subscript while holding the remaining parameters fixed at their true values. Let $\dot{l}_{\theta}$ be the ordinary score for $\theta$ when $S$ and $h$ are fixed. Then the efficient score function $\tilde{l}_{\theta}\in\{L_{2}^{0}(P_{0})\}^{p}$ for $\theta$ in the full model $\mathcal{P}$ at $P_{0}$ is $\tilde{l}_{\theta}=\dot{l}_{\theta}-\Pi_{0}(\dot{l}_{\theta}\mid\dot{\mathcal{P}}_{S}+\dot{\mathcal{P}}_{h})$, where $\Pi_{0}(l\mid\mathcal{W})$ denotes the orthogonal projection of $l$ onto the linear span of $\mathcal{W}$ [1].

2.1 Inference for Forward and Backward Recurrence Times

We now calculate the efficient score and information using the recurrence time $\tilde{T}$ potentially subject to right censoring, which covers both the forward and backward recurrence time settings. For the $i$-th individual, we observe ($\tilde{T}_{i}\wedge\tilde{C}_{i},\delta_{i},\tilde{Z}_{i}$), where $\tilde{T}_{i}$ is the recurrence time of the $i$-th individual, $\tilde{C}_{i}$ is the time of right censoring, which is assumed to be independent of the recurrence time conditionally given the covariates, $\delta_{i}$ is the indicator of whether the event time is observed and $\tilde{Z}_{i}$ is the $p\times 1$ observed covariate. Theorem 2.1 below demonstrates that in this setting, the efficient score equals that of the naive efficient estimator based on the conditional likelihood given the covariates.

For right-censored $\tilde{T}$, we assume $\Theta$ is a compact set in $\Re^{p}$, and that $\theta_{0}$ belongs in the interior of $\Theta$. For fixed but arbitrary $\theta$, we define our semiparametric model in terms of the distribution of $U(\theta)=e^{-\theta^{\prime}\tilde{Z}}\tilde{T}=e^{-(\theta-\theta_{0})^{\prime}\tilde{Z}}\tilde{U}$ and the corresponding censored variable $U^{c}(\theta)=e^{-\theta^{\prime}\tilde{Z}}\tilde{C}$. The conditional density of $U(\theta)$ given $\tilde{Z}=z$ is thus

[TABLE]

while the conditional hazard is

[TABLE]

Given $\tilde{Z}$ the density of $U(\theta)$ is monotone decreasing. We now state our assumptions:

• A1:

$\tilde{T}$ and $\tilde{C}$ are independent given $\tilde{Z}$;

• A2:

The distribution of $\tilde{C}$ is independent of the parameters $(\theta,S,h)$, and the distribution of $Z$ is independent of the parameters $(\theta,S)$;

• A3:

$\int S(v)dv<\infty$;

• A4:

$E_{g_{U(\theta)}}\left\{U^{2}\lambda(U)\right\}=\int u^{2}g_{U(\theta)}^{2}S(u)^{-1}du<\infty$.

The last assumption is needed to ensure that the density of $U(\theta)$ has finite Fisher information about $\theta$. The next theorem gives that the efficient score equals that from the naive efficient estimator.

Theorem 2.1**.**

Suppose that the covariate vector $\tilde{Z}$ is almost surely bounded. Define

[TABLE]

and

[TABLE]

Then under (A1)–(A4) and with $\phi(u)=1-ug(u)/S(u)$, the ordinary score for $\theta$ at $\theta=\theta_{0}$ is

[TABLE]

the tangent space $\dot{\mathcal{P}}_{S}$ for $S$ is $\{\dot{l}_{S}b:b\in L_{2}^{0}(S)\}$ where the score operator $\dot{l}_{S}$ for $S$ is given by

[TABLE]

the tangent space for $h$ is $\{k:k\in L_{2}(h),\int k(z)e^{\theta_{0}^{\prime}z}h(z)dz=0\}$, and the efficient score for $\theta$ at $\theta=\theta_{0}$ is

[TABLE]

Proof.

The likelihood for one observation $(U_{i}\wedge U^{c}_{i},\delta_{i},\tilde{Z}_{i})$ is given by

[TABLE]

Taking log and differentiating with respect to $\theta$, we obtain the ordinary score for $\theta$ at $\theta=\theta_{0}$,

[TABLE]

The expression in 2.3 for the ordinary score function for $\theta$ can be derived by noting that the quantity in brackets on the right hand side of the above expression is a stochastic integral with respect to the counting-process martingale in 2.2 [1], using proposition A.3.6 in [1].

Next, we can conclude from the Lemma 2.1 below that the tangent space $\dot{\mathcal{Q}}_{S}$ for $S$ can be considered the maximal tangent space $L_{2}^{0}(S)$. Hence the tangent space for $S$ can be expressed through the one dimensional submodels $\eta\mapsto S_{\eta}(t)=(1+\eta b(t))S(t)$ for any $b\in L_{2}^{0}(S)$, which yield the one dimensional baseline hazard submodels

[TABLE]

Differentiating the likelihood with respect to $\eta$ and setting $\theta=\theta_{0}$ now yields the score given in 2.4. In order to find $\Pi_{0}(\dot{l}_{\theta_{0}}\rvert\dot{\mathcal{P}}_{S})=\dot{l}_{S}b^{*}$ we find $b^{*}\in L^{0}_{2}(S)$ such that $\dot{l}_{\theta_{0}}-\dot{l}_{S}b^{*}\ \perp\ \dot{l}_{S}b$ for all $b\in L_{2}^{0}(S)$. That is $E\left\{\left(\dot{l}_{\theta_{0}}-\dot{l}_{S}b^{*}\right)\dot{l}_{S}b\right\}=0$. Note that $\dot{l}_{\theta_{0}}-\dot{l}_{S}b^{*}=\int_{-\infty}^{U^{c}(\theta_{0})}(\tilde{Z}R\phi-Rb^{*})dM(s)+(\tilde{Z}-E\tilde{Z})$. Conditioning on $\tilde{Z}$ and $U^{c}(\theta_{0})$ and using the fact that $U(\theta_{0})$ is distributed independently of $\tilde{Z}$ and $U^{c}(\theta_{0})$ we obtain

[TABLE]

The second equality above is obtained by using the result that if $Y_{i}=\int f_{i}dM,i=1,2$, then

[TABLE]

Thus $E\left\{\left(\dot{l}_{\theta_{0}}-\dot{l}_{S}b^{*}\right)\dot{l}_{S}b\right\}=0$ for all $b\in L_{2}^{0}(S)$ if

[TABLE]

Thus the projection of $\dot{l}_{\theta}$ on $\dot{\mathcal{P}}_{s}$ is given by

[TABLE]

Now for finding $\dot{\mathcal{P}}_{h}$ for $h\in\mathcal{H}^{\prime}$, we consider the one-parameter path $\eta\mapsto h_{\eta}=(1+\eta k)h$, where $k\in L_{2}^{0}(h)$. The score operator for $h$ is given by

[TABLE]

for $k\in L_{2}^{0}(h)$. Note that $\int k(z)e^{\theta_{0}^{\prime}z}h(z)dz=0$. If $h$ is unrestricted then the tangent space can be taken to be the orthocomplement of the linear span of $e^{\theta_{0}^{\prime}z}$, i.e., $[e^{\theta_{0}^{\prime}z}]^{\perp}$ in $L_{2}(h)$. Since $U(\theta_{0})$ is distributed independently of $\tilde{Z}$ and $U^{c}(\theta_{0})$, $E_{0}\{k(\tilde{Z})\dot{l}_{S}b\}=0$ for any $k\in[e^{\theta_{0}^{\prime}z}]^{\perp}$ and $\dot{l}_{S}b\in\dot{\mathcal{P}}_{S}$, i.e., $\dot{\mathcal{P}}_{S}\perp\dot{\mathcal{P}}_{h}$. Since $(z-E\tilde{Z})\in[e^{\theta_{0}^{\prime}z}]^{\perp}$, we obtain

[TABLE]

Now replacing $z$ with $\tilde{Z}$ and subtracting 2.6 and 2.7 from 2.3 yields the efficient score given in 2.5. The accelerated failure time model for $\tilde{T}$ given $\tilde{Z}$ is equivalent to the log-linear model $Y=\log(\tilde{T})=-\theta^{\prime}\tilde{Z}+\epsilon$, where $\epsilon$ has hazard function

[TABLE]

and $\lambda_{\tilde{T}}(u)$ is the baseline hazard for $\tilde{T}$. In the current setting,

[TABLE]

This model is the same as the linear regression model for $Y$ but with a sign change on $\theta$. The efficient score for the linear regression model under right-censoring is given in Expression (27) on Page 149 of [1] and has the same form as 2.5, except for changes in parameter and variable notation. Specifically, the function $R\phi(u)$ in 2.5 equals the negative of $R\psi(t)=-\dot{\lambda}(t)/\lambda(t)$ defined in Expression (23) of [1], after replacing $t$ with $\log(u)$, where the negative is due to the sign change. To see this, note that

[TABLE]

The first row follows from (2.8), the second row follows from the substitution $u=e^{t}$ followed by (2.9), and the last row follows from the definitions of $R$ and $\phi$. ∎

Thus the efficient score is free of $h$, so to estimate $\theta$ efficiently, one does not need to estimate the covariate distribution. Hence, one does not need to impose an additional identifiability condition for $h$ such as the mean-zero assumption. The efficient information is

[TABLE]

where $D(\tilde{Z},C,\theta_{0},s)=[\tilde{Z}-E\{\tilde{Z}|\ U^{c}(\theta_{0})\geq s\}]$. This is somewhat complicated to estimate, but the approach described in Remark 2 of [17] will yield a consistent estimator which can be used for inference on $n^{1/2}(\hat{\theta}_{n}-\theta_{0})$.

Since the backward recurrence times are uncensored, we can assume that the censoring times are infinite with $M(t)=I\{U(\theta)\leq t\}$ and $Ra(t)=a(t)$. Thus the efficient score for the backward recurrence time simplifies to

[TABLE]

and the efficient information becomes

[TABLE]

Before presenting Lemma 2.1, we provide a few needed definitions. Let ${\cal Q}$ be the model consisting of densities on $\mathcal{R}^{+}$ of the form $S_{g}(u)/\int_{0}^{\infty}S_{g}(v)dv$, where $g\in{\cal G}$; and let $\dot{\cal Q}_{S}$ and $\dot{\cal G}_{g}$ be the respective tangent sets for ${\cal Q}$ and ${\cal G}$ at $S$ and at $g$, where $g$ satisfies $S_{g}=S$. Lemma 2.1 establishes that $\dot{Q}_{S}=L_{2}^{0}(S)$, which is needed in the proof of Theorem 2.1 to identify $\dot{\mathcal{P}}_{S}$, a key technical step. We will be using score operators which allow us to construct scores for a model of interest from scores for a simpler model [see, e.g., Chapter 18 of 12].

Lemma 2.1**.**

If $A_{S}$ is the score operator mapping tangents in $\dot{\mathcal{G}}_{g}$ to $\dot{\mathcal{Q}}_{S}$, then $A_{S}\dot{\mathcal{G}}_{g}$ is dense in the maximal tangent set $L_{2}^{0}(S)$ for ${\cal Q}_{S}$, i.e., $\dot{\cal Q}_{S}=L_{2}^{0}(S)$.

Proof.

Let $g$ be the density on $\mathcal{R}^{+}$ corresponding to $S$. Consider the following parametric path through $g$:

[TABLE]

where $\psi:\mathcal{R}\mapsto\mathcal{R}^{+}$ is bounded, continuously differentiable with bounded derivative $\psi^{\prime}$ satisfying $\psi(0)=\psi^{\prime}(0)=1$ and $a\in L_{2}^{0}(g)$. Note that $L_{2}^{0}(g)$ is the closure within $L_{2}(g)$ of the derivatives of curves $g_{\eta}$ with respect to $\eta$ and $L_{2}^{0}(S)$ is the closure within $L_{2}(S)$ of derivatives of curves $\log(p_{g_{\eta}})$ with respect to $\eta$, where

[TABLE]

and $S_{\eta}$ is the survival function corresponding to $g_{\eta}$. Thus, $L_{2}^{0}(g)$ is the maximal non-parametric tangent set for $\mathcal{G}$ while $L_{2}^{0}(S)$ is the maximal tangent set for $\mathcal{Q}_{S}$. The corresponding parametric submodel for $p_{g}$ is

[TABLE]

Thus, the tangent set $\dot{\mathcal{Q}}_{S}$ (which consists of scores with respect to the one parameter models $p_{g_{\eta}}$) is given by the operator

[TABLE]

Let ${\cal A}$ be space of the bounded functions on $\mathcal{R}^{+}$ and ${\cal B}$ be the subset of ${\cal A}$ of functions which attain zero at all time points large enough. It is easy to verify that ${\cal A}$ is dense in $L_{2}^{0}(g)$ and that ${\cal B}$ is dense in $L_{2}^{0}(S)$, and, moreover, that $A_{S}a\in L_{2}^{0}(S)$ for all $a\in{\cal A}\cap L_{2}^{0}(g)$ and that $A_{S}^{\ast}b\in L_{2}^{0}(g)$ for all $b\in{\cal B}\cap L_{2}^{0}(S)$, where $A_{S}^{\ast}$ is the adjoint of $A_{S}$ defined as the solution to

[TABLE]

for all $a\in{\cal A}\cap L_{2}^{0}(g)$ and $b\in{\cal B}\cap L_{2}^{0}(S)$. This relation yields that $A_{S}^{\ast}b=\int_{0}^{u}b(v)dv$.

By definition of $A_{S}$, $\dot{\mathcal{Q}}_{S}$ is the closed linear span of $A_{S}{\cal A}$ in $L_{2}^{0}(S)$, and thus $\dot{\mathcal{Q}}_{S}\subset L_{2}^{0}(S)$. To prove the lemma, we need to verify that $L_{2}^{0}(S)\subset\dot{\mathcal{Q}}_{S}$ also holds. Suppose there is a $b_{0}\in L_{2}^{0}(S)$ which is not in $\dot{\mathcal{Q}}_{S}$. Then there exists a sequence $\{b_{n}\}\in{\cal B}\cap L_{2}^{0}(S)$ such that $\|b_{n}-b_{0}\|_{L_{2}^{0}(S)}\rightarrow 0$ and $0=$

[TABLE]

This now implies that $\int_{0}^{\infty}\left\{\int_{0}^{u}b_{n}(v)dv\right\}^{2}g(u)du\rightarrow 0$. We can now show that for any $0<c<\infty$ for which $S(c)>0$,

[TABLE]

as $n\rightarrow\infty$ by previous arguments combined with some analysis. Since $c$ was an arbitrary choice for which $S(c)>0$, we obtain that $b_{0}=0$ $S$-almost surely, and the desired conclusion follows. ∎

2.2 Inference for Length-Biased Data

A similar result may be obtained for length-biased data by replacing in the proof of Theorem 2.1 $S(u)/\int{S(u)d(u)}$ with $\tilde{g}=ug(u)/\int{S(u)d(u)}$, where, as before, $g$ is the density generating $S$. This yields the following result:

Theorem 2.2**.**

Using the same notation as Theorem 2.1 and under the same conditions, the efficient score for $\theta$ at $\theta=\theta_{0}$ for length-biased data is

[TABLE]

where for $a\in L_{2}^{0}(S)$,

[TABLE]

with $\phi(u)=1-ug(u)/S(u)$,

[TABLE]

and

[TABLE]

Proof.

The proof of Theorem 2.2 is very similar to the proof of Theorem 2.1. It follows along the same lines with a few minor differences, which are outlined below: The likelihood for one observation $(U_{i}\wedge U^{c}_{i},\delta_{i},\tilde{Z}_{i})$ is given by

[TABLE]

Here, the actual form of $g_{U(\theta)}(U_{i})=e^{(\theta-\theta_{0})^{\prime}\tilde{Z}_{i}}\tilde{g}(e^{(\theta-\theta_{0})^{\prime}\tilde{Z}_{i}}U_{i})$ is different from Theorem 2.1. Thus we need to replace ${\cal P}_{S}$ with ${\cal P}_{\tilde{g}}$, where the map $g\mapsto\tilde{g}$ is as implicitly defined above just before the statement of Theorem 2.2. Specifically, this is the new model for $g$ holding $\theta$ and $h$ fixed at their true values. The other models and submodels are the same as for Theorem 2.1 except that $S$ and $S_{0}$ are replaced by $\tilde{g}$ and $\tilde{g}_{0}$. Taking log of $l(\theta)$ and differentiating with respect to $\theta$ we obtain $\dot{l}_{\theta_{0}}=\tilde{Z}_{i}\left[\delta_{i}\phi(U_{i})+(1-\delta_{i})E\left\{\phi(U(\theta_{0}))\mid U(\theta_{0})>U^{c}_{i}\right\}\right]+(\tilde{Z}_{i}-E\tilde{Z})$ as the ordinary score for $\theta$ at $\theta=\theta_{0}$. The quantity in brackets on the right hand side is a stochastic integral with respect to the counting-process martingale in (2.14) and is thus also a martingale. Using this, we can obtain the ordinary score

[TABLE]

Let ${\cal G}$, $\dot{\cal G}_{g}$, and the model ${\cal Q}^{\prime}$ be as defined in Section 2 in the main text. Let $\dot{Q}_{\tilde{g}}^{\prime}$ be the tangent set for ${\cal Q}^{\prime}$ at $\tilde{g}$, and let $\dot{\cal P}_{\tilde{g}}$ be the tangent set corresponding to the new model ${\cal P}_{\tilde{g}}$. By using Lemma 2.2 below, which is similar to Lemma 2.1 but adapted to length-biased data, we can conclude that the tangent space $\dot{\mathcal{Q}}_{\tilde{g}}^{\prime}$ can be taken to be the maximal tangent space $L_{2}^{0}(\tilde{g})$, and thus we obtain that the score operator $\dot{l}_{\tilde{g}}$ for $\tilde{g}$ is

[TABLE]

In order to find $\Pi_{0}(\dot{l}_{\theta_{0}}\rvert\dot{\mathcal{P}}_{\tilde{g}})=\dot{l}_{\tilde{g}}b^{*}$, we find $b^{*}\in L^{0}_{2}(\tilde{g})$ such that $\dot{l}_{\theta_{0}}-\dot{l}_{\tilde{g}}b^{*}\ \perp\ \dot{l}_{\tilde{g}}b$ for all $b\in L_{2}^{0}(\tilde{g})$. That is $E\left\{\left(\dot{l}_{\theta_{0}}-\dot{l}_{\tilde{g}}b^{*}\right)\dot{l}_{\tilde{g}}b\right\}=0$. Note that $\dot{l}_{\theta_{0}}-\dot{l}_{\tilde{g}}b^{*}=\int_{-\infty}^{U^{c}(\theta_{0})}(\tilde{Z}R\phi-Rb^{*})dM(s)+(\tilde{Z}-E\tilde{Z})$.

After this, the proof in Theorem 2.1 can be followed verbatim to obtain the desired result. ∎

We now present Lemma 2.2 required for the proof of Theorem 2.2. This is a modification of Lemma 2.1 for the length-biased setting.

Lemma 2.2**.**

Consider the semi-parametric model $\mathcal{P}=\{P_{g}:g\in\mathcal{G}\}$, where the distribution $P_{g}$ has density $p_{g}(u)=ug(u)/\int S_{g}$ and $\mathcal{G}$ is a collection of densities on $\mathcal{R}^{+}$. Let $\dot{\mathcal{G}}_{g}$ and $\dot{\mathcal{P}}_{g}$ be the tangent sets for the models $\mathcal{G}$ and $\mathcal{P}$ respectively at $g$. If $A_{g}$ is the score operator mapping tangents in $\dot{\mathcal{G}}_{g}$ to $\dot{\mathcal{P}}_{g}$ then, $A_{g}\dot{\mathcal{G}}_{g}$ is dense in the maximal tangent set $L_{2}^{0}(S)$ for $\mathcal{P}$.

Proof.

Consider the following parametric path through $g$:

[TABLE]

where $\psi:\mathcal{R}\mapsto\mathcal{R}^{+}$ is bounded, continuously differentiable with bounded derivative $\psi^{\prime}$ satisfying $\psi(0)=0$ and $\psi^{\prime}(0)=g(0)$ and $a\in L_{2}^{0}(g)$. Thus, using similar notations as in Lemma 1, $L_{2}^{0}(g)$ is the maximal non-parametric tangent set for $\mathcal{G}$ while $L_{2}^{0}(\tilde{g})$ is the maximal tangent set for $\mathcal{Q}^{\prime}$. The corresponding parametric submodel for $\tilde{g}$ will then be

[TABLE]

After this, the proof of Lemma 2.1 can be followed essentially verbatim to obtain the desired result. ∎

3 Simulation Studies

As our method does not rely on the covariate distribution, it is a special case of the model used in the paper by [17] (since we assume that the covariates are constant over time). So, we use their profile likelihood approach to estimate $\theta$ and compare it with Klaassen’s mean zero approach and also the known covariate structure approach. We consider only 1 covariate Z $\sim\mbox{Unif}(-1,1)$. So, $\tilde{Z}$ has density given by $\theta e^{\theta z}/(e^{\theta}-e^{-\theta})$, where $-1\leq z\leq 1$. We take different values of $\theta$ and assume that the error distribution is standard normal, i.e., $U$ is lognormal. Then, we use all three methods to estimate $\theta$. For the profile-likelihood approach, we use the Gaussian kernel and a bandwidth of $h_{n}=Qn^{-1/5}$, where Q is the interquartile range of the data. We consider 1000 replicates and look at the mean bias and variance in estimating $\theta$. We also look at what happens when the covariate distribution is misspecified. For this, we consider Z$\sim\mbox{Unif}(x,1)$ for some choice of x. So, $\tilde{Z}$ has density given by $\theta e^{\theta z}/(e^{\theta}-e^{-\theta x})$ where x$\leq$z$\leq$1. We take the values $x=-0.9$ and $-0.8$ and compare both the mean zero and known covariate distribution, assuming Z $\sim\mbox{Unif}(-1,1)$. We take $\theta$=1 for these simulations. We consider 1000 replicates in this case as well. The results are given in Table 1.

Thus we find that the estimates obtained using our methods are quite comparable to the special case where the covariance structure is known, although Klaassen’s method has lower variance. This is expected because Klaassen’s method is under some additional model assumptions which are not used by our method. However, their estimates are very sensitive to model specification. On the other hand, our naive analysis yields unbiased estimates in both cases. The variance estimators accurately reflect the actual variance, while the confidence intervals also have correct coverage probabilities.

4 Data Analysis

For illustration, we analyze a subset of the data used by [11]. It is a backward recurrence time data setting on the time to pregnancy obtained from a large French telephone survey. Women were eligible if they were between 18-44 years old, were living with a male partner and did not use any method to avoid pregnancy. We consider only nulliparous women who had not initiated any fertility treatment. The response variable was the current duration of unprotected intercourse, which is the time elapsed from the start of unprotected intercourse and the interview. The estimates obtained for the covariates along with the 95 $\%$ confidence intervals are given in Table 2. We note that the naive estimator can accurately determine the effect of the covariates and is comparable with the ordinary least squares results.

5 Discussion

The assumption that the initiating event follows a homogeneous Poisson process has been widely used [6, 14, 9]. Such an assumption may be reasonable over short time periods, where the rate of the initiating event may be stable. Over longer time periods, where the assumption may not be valid, alternative methods may be needed. This is a challenging problem worthy of further investigation.

In the presence of covariates, a popular alternative to the accelerated failure time model is the proportional hazards model [7] given by $\lambda_{T\mid Z}(t)=e^{\theta^{\prime}Z}\lambda(t)$, where $\lambda_{T\mid Z}$ is the hazard function of $T$ given the covariate vector $Z$ and $\lambda$ is an unspecified baseline hazard function. Here the density of $T$ is given by $e^{\theta^{\prime}z}\lambda(t)e^{-e^{\theta^{\prime}z}\Lambda(t)}$, where, $\Lambda$ is the cumulative baseline hazard function satisfying $\Lambda(0)=0$. If we assume a proportional hazards model for $T$, then by (1.1), under length-biased and cross-sectional sampling, the conditional density of the forward or the backward recurrence time $\tilde{T}$, given $\tilde{Z}$, is

[TABLE]

and the conditional hazard function is

[TABLE]

Thus, going from $T$ to $\tilde{T}$, the proportional hazard structure is lost, unless either the baseline hazard is constant or $T$ given $Z$ follows a Pareto distribution [14]. The usual techniques for fitting the proportional hazards model may yield biased results. Furthermore, it is unclear whether information in the observed covariates may be employed to yield more efficient estimation, in contrast to our results for the accelerated lifetime model for T given Z. Analogous issues arise for the length-biased time $T_{\mathrm{LB}}$ under the proportional hazards model. Further investigation is needed in these areas.

Under additional assumptions on the covariate distribution, for example, $Z$ is mean zero or $h$ is known or specified parametrically, information may be gained using the observed covariates [14]. The trade-off for these efficiency gains is the potential bias associated with the additional modeling assumptions for the distribution of $Z$. In our simulations given above, we found that the efficient estimators are rather sensitive to the extra assumptions and that the gains in efficiency when such assumptions are correctly specified are rather modest compared to the naive efficient estimators.