Censored Interquantile Regression Model with Time-Dependent Covariates

Abstract

Conventionally, censored quantile regression stipulates a specific, pointwise conditional quantile of the survival time given covariates. Despite its model flexibility and straightforward interpretation, the pointwise formulation oftentimes yields rather unstable estimates across neighboring quantile levels with large variances. In view of this phenomenon, we propose a new class of quantile-based regression models with time-dependent covariates for censored data. The models proposed aim to capture the relationship between the failure time and the covariate processes of a target population that falls within a specific quantile bracket. The pooling of information within a homogeneous neighborhood facilitates more efficient estimates hence, more consistent conclusion on statistical significances of the variables concerned. This new formulation can also be regarded as a generalization of the accelerated failure time model for survival data in the sense that it relaxes the assumption of global homogeneity for the error at all quantile levels. By introducing a class of weighted rank-based estimation procedure, our framework allows a quantile-based inference on the covariate effect with a less restrictive set of assumptions. Numerical studies demonstrate that the proposed estimator outperforms existing alternatives under various settings in terms of smaller empirical biases and standard deviations. A perturbation-based resampling method is also developed to reconcile the asymptotic distribution of the parameter estimates. Finally, consistency and weak convergence of the proposed estimator are established via empirical process theory. Supplementary materials for this article are available online.

KEYWORDS:

• Accelerated failure time model
• Censored data
• Empirical process
• Quantile inference
• Survival analysis
• Time-dependent covariate

1 Introduction

Censored quantile regression has been a popular alternative to Cox proportional hazards model (Cox 1972, 1975) in survival analysis. This broad class of models allows a greater degree of flexibility as it can capture the heterogeneity of covariate effects without assuming that the relative risks are constant with time. It also provides better interpretability as the model directly relates a specific, pointwise quantile of the failure time with a set of candidate covariates. There has been rich literature discussing quantile regression models on censored survival data. In particular, Powell (1984, 1986) studied the least absolute deviation estimation for the censored regression model. Under the independent censoring assumption, Ying, Jung, and Wei (1995) derived the median regression procedure for survival data. With a less restrictive assumption on conditionally independent censoring, Portnoy’s (2003) and Peng and Huang’s (2008) proposals provided novel adaptations for Kaplan-Meier (KM) and Nelson-Aalen (NA) estimators, respectively, upon which corresponding inference procedures for censored quantiles can be implemented. To avoid recursive algorithms required in the aforementioned procedures, Wang and Wang (2009) insightfully developed a locally weighted procedure that adopts the redistribution-of-mass idea for a kernel-based, local reweighting scheme. All these approaches were, however, designed for time-independent covariates only.

In many real applications, since survival events occur over time, there is a natural need to handle time-dependent covariates such as heart blood pressures or the drug dosages received by individual patients. Statistical theories of the Cox proportional hazards model can be readily extended to handle time-dependent covariates, see, for example, Andersen and Gill (1982), Kalbfleisch and Prentice (2002) and Martinussen and Scheike (2006) amongst others. The classic accelerated failure time model has also been extended in Lin and Ying (1995), Zeng and Lin (2007), and Tseng, Hsieh, and Wang (2005), amongst others to accommodate time-dependent covariates. For censored quantile regression with time-dependent covariates, there have only been few results to the best of our knowledge. One interesting work is Gorfine, Goldberg, and Ritov (2017), whose estimating equation was developed upon an inverse probability weighting scheme under the independent censoring assumption.

As pointed out in Peng and Huang (2008), censored quantile regression models can be reduced to the accelerated failure time model (Robins and Tsiatis 1992; Lin and Ying 1995) if the regression parameters are constant for all quantile values $\tau \in (0,1)$. In other words, when the homogeneity assumption holds for all τ’s, the two models are equivalent. This observation motivates us to explore the uncharted region between the two ends of the spectrum. Indeed, our new proposal attempts to seek a reasonable compromise between model flexibility and estimation efficiency. Because of the homogeneity assumption, the accelerated failure time model can pool the information provided by all the observations for estimating a reduced set of model parameters. As we can see in Example 4 in Section 5, if one can delicately leverage the local homogeneity shared amongst the specific quantile level concerned and its neighboring levels, the problem of obtaining volatile point estimates for quantile regression parameters, which is frequently observed in many existing studies, can be alleviated.

Recall that the censored quantile regression model studies the influence of the covariates on the failure time at designated quantile levels, such as the median and the quartiles, while the accelerated failure time model examines the overall covariate effect provided the aforementioned homogeneity assumption holds for all quantile levels. None of these models, however, can quantify the covariate effects over a specific group in the sample, say the top 25% or the bottom 25% of a population in term of length of survival after covariate adjustment. For instance, to measure the effectiveness of a new drug, investigators are interested in quantifying its effect on the survival times of individuals who are among the shortest lived 25% in the population given the covariates instead of its effect on an individual who represents the lower quartile of the failure time distribution. In such scenarios, a localized censored quantile regression model with parameters representing the treatment effect over the designated range of quantiles should appear more attractive, whereas the classic quantile regression model only quantifies treatment effect on particular levels of quantiles instead.

The main contribution of this work is the development of a quantile-based version of the accelerated failure time model for censored data with time-dependent covariates under which we relax the homogeneity requirement under the conditional independent censoring assumption. Moreover, our approach also provides an alternative perspective to the study of quantile regression problems. The pooling effect enjoyed by our model can remediate the drawback that many existing proposals for censored quantile regression may share. Besides numerical stability, our model also offers an alternative perspective for researchers so as to best match their primary goals of quantifying the covariate effects on the failure time for an interval of quantiles rather than a specific level. It should be emphasized that although there are similar works in literature that discuss methods of combining neighboring quantile coefficient estimates to yield more stable estimates, to the best of our knowledge, none can be straightforwardly extended to our specific problem. In Zheng, Peng, and He (2015), for example, the aggregation is carried out via the penalty term based on an overall sum of the magnitudes of the neighboring coefficients concerned, but the parameter estimates remain capturing the covariates effects on the response variable at specific quantile levels. Same as Zheng, Peng, and He (2015), Zou and Yuan (2008), and Kai, Li, and Zou (2010) are also designed for uncensored observations in which case the celebrated check function can be readily applied to provide initial estimates upon which smoothing or other aggregation is introduced to achieve the inference goals.

The rest of this article is organized as follows: Section 2 outlines the model framework and establishes the corresponding estimation equation for the model parameters. Section 3 outlines the asymptotic results of the proposed estimator, while Section 4 provides a perturbation based resampling method and a two stage testing procedure for inference. A set of simulation results is presented in Section 5, followed by a data analysis of the Stanford heart transplant data in Section 6. A short concluding discussion is given in Section 7. All proofs of the asymptotic results and justification of the resampling method are presented in the appendix. A sample of the MATLAB code used in our numerical studies is provided online as a supplementary material.

2 Methodology

For $i=1,\dots ,n$, let T_i denote the actual failure time of the ith individual, ${\overline{\mathbf{Z}}}_i=\{{\mathbf{Z}}_i(s),s\ge 0\}=\{{(Z_{i1}(s),\dots ,Z_{ip}(s))}^{\top },s\ge 0\}$ be the p-dimensional multivariate random covariate process corresponding to this person, and ${\mathit{\beta }}_0$ be the p-dimensional vector of coefficients. Let ${\epsilon }_i=h_i(T_i;{\mathit{\beta }}_0)$ be a monotone transform of T_i, where $$h_i(T_i;\mathit{\beta })={\int }_0^{T_i}\hspace{0.17em}\exp \hspace{0.17em}\{-{\mathbf{Z}}_i{(s)}^{\top }\mathit{\beta }\} ds.$$

Note that $h_i(T_i;{\mathit{\beta }}_0)$ reduces to $T_i\hspace{0.17em}\exp \hspace{0.17em}\{-{\mathbf{Z}}_i^{\top }{\mathit{\beta }}_0\}$ when ${\mathbf{Z}}_i$ is time-independent, that is, we have $\hspace{0.17em}\log \hspace{0.17em}{\epsilon }_i=\hspace{0.17em}\log \hspace{0.17em}{T}_i-{\mathbf{Z}}_i^{\top }{\mathit{\beta }}_0$, which reduces to the accelerated failure time model with time-independent covariates. In this regard, the definition of ${\epsilon }_i$ here generalizes the error term to the setting of time-dependent covariates, which can be interpreted as a rescaled failure time of the ith individual such that these transformed time stamps are identically distributed under the accelerated failure time model.

In our setup, instead of assuming a homogeneous error term, that is, a global homogeneity for all quantile levels, we hypothesize only a local homogeneity for a designated range of quantiles. More precisely, our model assumes that for all $u\in [{\tau }_L,{\tau }_U]\subset [0,1]$, the uth quantile of ${\epsilon }_i$ conditional on the covariate process ${\overline{\mathbf{Z}}}_i$ is homogeneous for all $i=1,\dots ,n$. Mathematically, it stipulates that (2.1) $$Q_{h_i(T_i;{\mathit{\beta }}_0)}(u\mid {\overline{\mathbf{Z}}}_i)=q(u),\hspace{1em}(u\in [{\tau }_L,{\tau }_U];\hspace{1em}i=1,\dots ,n),$$(2.1) where q(u) is a function independent of ${\overline{\mathbf{Z}}}_i$ representing the common conditional quantile, and $Q_{h_i(T_i;{\mathit{\beta }}_0)}(·\mid {\overline{\mathbf{Z}}}_i)$ denotes the conditional quantile function of $h_i(T_i;{\mathit{\beta }}_0)$ given ${\overline{\mathbf{Z}}}_i$. Our goal is to estimate ${\mathit{\beta }}_0\equiv {\mathit{\beta }}_{0,[{\tau }_L,{\tau }_U]}$ under model (2.1), where ${\mathit{\beta }}_0$ quantifies the covariate effect over the homogeneous region. Here, we suppress the dependence on τ_L and τ_U for brevity of notation. Note that the choice of $[{\tau }_L,{\tau }_U]$ depends on the primarily interested range of quantiles. As an example, if the user is interested in the covariate effects on the bottom 25% subgroup of population after covariate adjustment, then $[{\tau }_L,{\tau }_U]$ could be chosen as $[0,0.25]$.

Model (2.1) can be viewed as a generalization of the accelerated failure time model with time-dependent covariates (Robins and Tsiatis 1992; Lin and Ying 1995). Recall their model assumes that $h_i(T_i;{\mathit{\beta }}_0)$ shares a common baseline failure time distribution conditional on ${\overline{\mathbf{Z}}}_i$, which is equivalent to say $Q_{h_i(T_i;{\mathit{\beta }}_0)}(u\mid {\overline{\mathbf{Z}}}_i)=q(u)$ for all $u\in [0,1]$, that is, (2.1) holds for ${\tau }_L=0$ and ${\tau }_U=1$. In contrast, we relax such a requirement so that (2.1) is only satisfied for a sub-interval of quantiles $[{\tau }_L,{\tau }_U]$. Recall that the accelerated failure time model can be regarded as a special case of the censored quantile regression model which does not assume any error homogeneity. In this regard, our model can be viewed as a hybrid that connects these two extremes. While the parameter in the accelerated failure time model gives the right transformation of the failure time so that the rescaled time stamps $h_i(T_i;{\mathit{\beta }}_0)$ are homogeneous, our model parameter can be regarded as an analogy that gives the correct transformation over $[{\tau }_L,{\tau }_U]$. In other words, one can interpret ${\mathit{\beta }}_0$ as the uniform covariate effect on the failure time over this subgroup of population that lies within this quantile bracket. As an example, suppose the homogeneity within the middle 50% of quantile levels is concerned, that is, $[{\tau }_L,{\tau }_U]=[0.25,0.75]$, then our model hypothesizes that a unit uniform increment in the covariate process ${\mathbf{Z}}_i(s)$ results in $\hspace{0.17em}\exp \hspace{0.17em}\left\{{\mathbf{1}}^{\top }{\mathit{\beta }}_0\right\}$-times prolongation of the failure time for this middle group of sub-population, where $\mathbf{1}={(1,\dots ,1)}^{\top }$.

Suppose our data are iid copies of $\{{\tilde{T}}_i,{\Delta }_i,{\overline{\mathbf{Z}}}_i({\tilde{T}}_i)\}$ for $i=1,\dots ,n$, where ${\tilde{T}}_i=T_i\wedge C_i$ is the observed event time, ${\Delta }_i=I(T_i\le C_i)$ is the censoring indicator, and ${\overline{\mathbf{Z}}}_i({\tilde{T}}_i)=\{{\mathbf{Z}}_i(s),0\le s\le {\tilde{T}}_i\}$ is the covariate history up to the survival time. Here, C_i is the censoring time which is conditionally independent of T_i given ${\overline{\mathbf{Z}}}_i$, and $I(·)$ denotes the indicator function. Define the conditional cumulative hazard function, counting process and at-risk process associated with $h_i(T_i;{\mathit{\beta }}_0)$ to be ${\Lambda }_i(t;{\mathit{\beta }}_0)$, $N_i(t;{\mathit{\beta }}_0)={\Delta }_{i}I(h_i({\tilde{T}}_i;{\mathit{\beta }}_0)\le t)$, $Y_i(t;{\mathit{\beta }}_0)=I(h_i({\tilde{T}}_i;{\mathit{\beta }}_0)\ge t)$, respectively. Note that model (2.1) implies that $h_i(T_i;{\mathit{\beta }}_0)$ has a homogeneous conditional distribution function for $t\in [q({\tau }_L),q({\tau }_U)]$, so we can write ${\Lambda }_i(t;{\mathit{\beta }}_0)=\Lambda (t;{\mathit{\beta }}_0)$ for such t. Now, based on the classical counting process theory (Fleming and Harrington 1991), we have (2.2) $$E\left[d{N}_i(t;{\mathit{\beta }}_0)-Y_i(t;{\mathit{\beta }}_0)d\Lambda (t;{\mathit{\beta }}_0)\mid {\overline{\mathbf{Z}}}_i\right]=0,$$(2.2) and also (2.3) $$E\left[{\mathbf{Z}}_i(h_i^{-1}(t;{\mathit{\beta }}_0))\{d{N}_i(t;{\mathit{\beta }}_0)-Y_i(t;{\mathit{\beta }}_0)d\Lambda (t;{\mathit{\beta }}_0)\}\mid {\overline{\mathbf{Z}}}_i\right]=0,$$(2.3) for $t\in [q({\tau }_L),q({\tau }_U)]$, where $h_i^{-1}(·;;\mathit{\beta })$ is the inverse function of $h_i(·;;\mathit{\beta })$. By solving the empirical counterpart of (2.2), we obtain a Nelson-Aalen type estimator of $d\Lambda (t;{\mathit{\beta }}_0)$, which is ${\sum }_{i=1}^{n}d{N}_i(t;{\mathit{\beta }}_0)/{\sum }_{i=1}^n{Y}_i(t;{\mathit{\beta }}_0)$ for $t\in [q({\tau }_L),q({\tau }_U)]$. Meanwhile, with a slight abuse of notation, we denote the τth unconditional quantile of $h_i(T_i;\mathit{\beta })$ by $q(\tau ;\mathit{\beta })$, that is, $q(\tau ;\mathit{\beta })$ is the inverse function of $F(t;\mathit{\beta })$, where $F(t;\mathit{\beta })$ is the unconditional, or marginal, distribution function of $h_i(T_i;\mathit{\beta })$. Model (2.1) implies the homogeneous conditional quantile $q(\tau )$ coincides with $q(\tau ;{\mathit{\beta }}_0)$ for $\tau \in [{\tau }_L,{\tau }_U]$, so for $\tau ={\tau }_L$ and $\tau ={\tau }_U,q(\tau )$ can be estimated by $\widehat{q}(\tau ;{\mathit{\beta }}_0)$, the inverse function of $\widehat{F}(t;{\mathit{\beta }}_0)$, where $\widehat{F}(t;\mathit{\beta })$ is an estimator of $F(t;\mathit{\beta })$. In particular, we take $\widehat{F}(t;\mathit{\beta })$ to be the Kaplan-Meier estimator based on $\{h_i({\tilde{T}}_i;{\mathit{\beta }}_0),{\Delta }_i\}$. As mentioned, this estimation works because the conditional quantile and the unconditional quantile are equivalent for $\tau \in [{\tau }_L,{\tau }_U]$ when $\mathit{\beta }={\mathit{\beta }}_0$.

Incorporate the two estimators of $d\Lambda (t;{\mathit{\beta }}_0)$ and $q(\tau )$ into the empirical counterpart of (2.3) and integrate over the range of t, it is natural to consider the weighted estimating equation $U_n(\mathit{\beta };{\tau }_L,{\tau }_U)=0$, where $$U_n(\mathit{\beta };{\tau }_L,{\tau }_U)=\frac{1}{n}\sum _{i=1}^n{\int }_{\widehat{q}({\tau }_L;\mathit{\beta })}^{\widehat{q}({\tau }_U;\mathit{\beta })}\Psi (t,\mathit{\beta })$$ $$\times \left\{{\mathbf{Z}}_i(h_i^{-1}(t;\mathit{\beta }))-\frac{S^{(1)}(t;\mathit{\beta })}{S^{(0)}(t;\mathit{\beta })}\right\}d{N}_i(t;\mathit{\beta }).$$

Here, $S^{(0)}(t;\mathit{\beta })=n^{-1}{\sum }_{i=1}^n{Y}_i(t;\mathit{\beta }),S^{(1)}(t;\mathit{\beta })=n^{-1}{\sum }_{i=1}^n$ ${\mathbf{Z}}_i(h_i^{-1}(t;\mathit{\beta }))Y_i(t;\mathit{\beta })$, and $\Psi (t,\mathit{\beta })$ is a weight function that possibly depends on the data $\{{\tilde{T}}_i,{\Delta }_i,{\overline{\mathbf{Z}}}_i({\tilde{T}}_i),i=1,\dots ,n\}$. Our proposed estimator of ${\mathit{\beta }}_0$, denoted by $\widehat{\mathit{\beta }}$, is obtained as the root of $U_n(\mathit{\beta };{\tau }_L,{\tau }_U)=0$. To this end, we simplify $U_n(\mathit{\beta };{\tau }_L,{\tau }_U)$ into a more familiar form by recalling that $N_i(t;\mathit{\beta })={\Delta }_{i}I(h_i({\tilde{T}}_i;\mathit{\beta })\le t)$. The above display can thus be rewritten as (2.4) $$\begin{array}{c}{U}_n(\mathit{\beta };{\tau }_L,{\tau }_U)=\frac{1}{n}\sum _{i=1}^n\Psi (h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })\\ \times \left\{{\mathbf{Z}}_i({\tilde{T}}_i)-\frac{S^{(1)}(h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })}{S^{(0)}(h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })}\right\}\\ \hspace{1em}\times \{N_i(\widehat{q}({\tau }_U;\mathit{\beta });\mathit{\beta })-N_i(\widehat{q}({\tau }_L;\mathit{\beta });\mathit{\beta })\}\mathrm{}.\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\end{array}$$(2.4)

Here, $\widehat{q}({\tau }_L;\mathit{\beta })$ reduces to 0 if ${\tau }_L=0$ and $\widehat{q}({\tau }_U;\mathit{\beta })$ reduces to $\infty$ if ${\tau }_U=1$. Observe that if $h_i({\tilde{T}}_i;\mathit{\beta })>h_j({\tilde{T}}_j;\mathit{\beta })$, then $Y_j(h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })=0$; otherwise if $h_i({\tilde{T}}_i;\mathit{\beta })\le h_j({\tilde{T}}_j;\mathit{\beta })$, then we have $h_j^{-1}(h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })\le {\tilde{T}}_j$. Therefore, both $S^{(1)}(h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })$ and $S^{(0)}(h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })$ are observable for any $\mathit{\beta }$ given the data, regardless of the censoring status of the subject. As a result, despite we condition on the entire covariate process ${\overline{\mathbf{Z}}}_i$ in the derivation, the estimating (2.4) remains computable provided that only ${\overline{\mathbf{Z}}}_i({\tilde{T}}_i)$ is known, that is, the covariate value up to the observed survival time instead of the entire process. We emphasize that, as in the accelerated failure time model and the quantile regression model for time-dependent covariates such as Lin and Ying (1995) and Gorfine, Goldberg, and Ritov (2017), although internal covariates could possibly be fit into (2.4), one has to be cautious in interpreting the model parameters as opposed to external covariates; see chap. 6 of Kalbfleisch and Prentice (2002).

Given that $U_n(\mathit{\beta };{\tau }_L,{\tau }_U)$ is neither continuous nor monotone in general, we solve its root by minimizing the L₁ norm. In particular, we adopt the differential evolution algorithm of Rainer Storn for root solving. The corresponding MATLAB code devec3.m is available in the author’s website at http://www.icsi.berkeley.edu/∼storn/code.html. According to our numerical experience, such stochastic algorithm is generally more robust to local minima than, for example, the deterministic simplex search method. Nevertheless, when the covariate is high-dimensional, the computational burden will increase substantially due to the curse of dimensionality. It is interesting for explore the corresponding variable selection issue for this problem, as in, for example, Zheng, Peng, and He (2018) in the future.

To shed light on the connection between our formulation with existing estimators, one could view $\Psi (h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })$ $\{N_i(\widehat{q}({\tau }_U;\mathit{\beta });\mathit{\beta })-N_i(\widehat{q}({\tau }_L;\mathit{\beta });\mathit{\beta })\}$ as a single subject specific weight, then $U_n(\mathit{\beta };{\tau }_L,{\tau }_U)$ falls into the class of weighted log rank statistics for time-dependent covariate setting, with $\Psi (h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })\{N_i(\widehat{q}({\tau }_U;\mathit{\beta });\mathit{\beta })-N_i(\widehat{q}({\tau }_L;\mathit{\beta });\mathit{\beta })\}$ as the modified weight. Indeed, when $\Psi (h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })=1$, ${\tau }_L=0$ and ${\tau }_U=1$, the expression reduces to Δ_i and thus (2.4) recovers the estimating equation of Lin and Ying (1995).

Typical weight functions include $\Psi (t;\mathit{\beta })=1$ and $\Psi (t;\mathit{\beta })=S^{(0)}(t;\mathit{\beta })$, which is generally known as Gehan weight. Due to the extra complication introduced by the quantile specific inference and the time-dependent covariates, the two choices do not allow for an easier numerical implementation as in Jin et al. (2003). Based on our numerical experience as shown in Section 5, the performance using the two different weight functions does not differ considerably. Unless otherwise specified, we shall assume the estimating equation is unweighted when we adopt (2.4) in the numerical studies.

Before concluding this section, we stress here that despite model (2.1) resembles the formulation considered in Gorfine, Goldberg, and Ritov (2017), they are fundamentally different in the following two senses. First, Gorfine, Goldberg, and Ritov (2017) developed their inference procedure based on inverse probability weighting approach because the event time and the censoring time concerned are assumed to be independent. In particular, they made use of the observation that $P(h_i(T_i;{\mathit{\beta }}_{\tau })\le 1\mid {\overline{\mathbf{Z}}}_i)=\tau$ to derive their estimating equation, namely $$\frac{1}{n}\sum _{i=1}^n\frac{{\Delta }_i{\mathbf{Z}}_i(0)}{\widehat{G}({\tilde{T}}_i)}\left\{I\left(h_i(T_i;{\mathit{\beta }}_{\tau })\le 1\right)-\tau \right\}=0,$$where $\widehat{G}(·)$ denotes the Kaplan-Meier estimate of the survival function of the censoring distribution. On the contrary, in this article, we relax this assumption and adopt the martingale formulation as in (2.2) through which we establish our inference procedure. Second, in Gorfine, Goldberg, and Ritov (2017), specific quantile parameters are estimated independently with no pooling of information across adjoining quantile levels. As we witness in our numerical studies, overlooking such homogeneity results in volatile estimates, see Example 4 in Section 5. This observation indeed motivates the development of our methodology.

Another important difference between our formulation and Gorfine, Goldberg, and Ritov (2017) is the exclusion of intercept term in the covariates. Note that their intercept parameter represents the quantile of the error term, which is identifiable because of their additional model assumption that requires this quantity to be one; see eq. (2.4)(2.4) $$\begin{array}{c}{U}_n(\mathit{\beta };{\tau }_L,{\tau }_U)=\frac{1}{n}\sum _{i=1}^n\Psi (h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })\\ \times \left\{{\mathbf{Z}}_i({\tilde{T}}_i)-\frac{S^{(1)}(h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })}{S^{(0)}(h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })}\right\}\\ \hspace{1em}\times \{N_i(\widehat{q}({\tau }_U;\mathit{\beta });\mathit{\beta })-N_i(\widehat{q}({\tau }_L;\mathit{\beta });\mathit{\beta })\}\mathrm{}.\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\end{array}$$(2.4) of Gorfine, Goldberg, and Ritov (2017). In our setup, the distribution of the error term is left to be entirely unspecified without any assumption on its quantile. In other words, we combine the intercept and the error into one quantity. Hence, we skip the intercept estimate in our discussion as our main focus remains primarily on the covariate effects.

3 Asymptotic Results

In this section, we state the asymptotic results of the proposed estimator and the corresponding regularity conditions, while the proofs are deferred to the appendix. The required conditions are listed as follows:

• Condition 1. $|Z_{ik}(0)|+{\int }_0^{\infty }|d{Z}_{ik}(t)|<\infty$ for each $i=1,\dots ,n$ and $k=1,\dots ,p$.

• Condition 2. ${\mathit{\beta }}_0\in \mathcal{B}$, where $\mathcal{B}$ is a bounded convex region.

• Condition 3. $P(\tilde{T}>L\mid \overline{\mathbf{Z}})>0$ almost surely for all $\overline{\mathbf{Z}}$, where $L=\underset{\mathit{\beta }\in \mathcal{B}}{\sup }{h}^{-1}(t;\mathit{\beta })$.

• Condition 4. $\underset{t\ge 0,\mathit{\beta }\in \mathcal{B}}{\inf }f(t;\mathit{\beta })>0$, where $f(t;\mathit{\beta })$ denotes the density function of the unconditional distribution of $h(T;\mathit{\beta })$.

• Condition 5. $\underset{\mathit{\beta }\in \mathcal{B}}{\sup }|\Psi (0;\mathit{\beta })|+{\int }_0^{\infty }|d\Psi (t;\mathit{\beta })|<\infty$, and there exists a function $\psi (t;\mathit{\beta })$ such that $\underset{t\in [q({\tau }_L;\mathit{\beta }),q({\tau }_U;\mathit{\beta })],\mathit{\beta }\in \mathcal{B}}{\sup }$ $|\Psi (t;\mathit{\beta })-\psi (t;\mathit{\beta })|\to 0$ almost surely.

• Condition 6. ${\mathbf{A}}_n({\mathit{\beta }}_0;{\tau }_L,{\tau }_U)$ is invertible, where ${\mathbf{A}}_n(\mathit{\beta };{\tau }_L,{\tau }_U)$ is the gradient of ${\tilde{U}}_n(\mathit{\beta };{\tau }_L,{\tau }_U)$ with respect to $\mathit{\beta }$, and ${\tilde{U}}_n(\mathit{\beta };{\tau }_L,{\tau }_U)$ is a smoothed analogy of $U_n(\mathit{\beta };{\tau }_L;{\tau }_U)$ as defined in the appendix.

Condition 1 is a mild smoothness requirement on the covariate processes. In particular, if $Z_{ik}(t)$ has a uniformly bounded derivative, then it is absolutely continuous and hence of bounded total variation. More generally, $Z_{ik}(t)$ is allowed to have jumps with sojourn times having uniformly bounded densities. Conditions 2, 3, and 4 are standard assumptions in survival analysis. In particular, using the idea of the functional delta method, see, for example, sec. 12 of Kosorok (2008), Conditions 3 and 4 imply the uniform consistency of the Kaplan-Meier estimator $\widehat{F}(t;\mathit{\beta })$ and hence the quantile estimator $\widehat{q}(\tau ;\mathit{\beta })$ over $\mathit{\beta }\in \mathcal{B}$. Condition 5 imposes restriction on the weight function so as to guarantee that it is well behaved to achieve the consistency and weak convergence of the proposed estimator. This condition is satisfied for a broad class of weight functions including, but not limited to, the constant and Gehan weights discussed previously. Condition 6 is essential to ensure the uniqueness of ${\mathit{\beta }}_0$ and to recover the asymptotic normality of the estimator from that of the estimating equation.

With the above conditions, we now establish the consistency and weak convergence of the proposed estimator through the following theorems.

Theorem 1.

If Conditions 1–6 hold, then there exists a fixed neighborhood of ${\mathit{\beta }}_0$, denoted by $\mathcal{N}({\mathit{\beta }}_0)$, such that for any $\widehat{\mathit{\beta }}\in \mathcal{N}({\mathit{\beta }}_0)$, we have $\Vert \widehat{\mathit{\beta }}-{\mathit{\beta }}_0\Vert \to 0$ almost surely.

Theorem 2.

If Conditions 1–6 hold, then $n^{1/2}(\widehat{\mathit{\beta }}-{\mathit{\beta }}_0)$ converges weakly to a mean zero multivariate normal with covariance matrix $\mathbf{A}{({\mathit{\beta }}_0;{\tau }_L,{\tau }_U)}^{-1}\mathbf{V}({\mathit{\beta }}_0;{\tau }_L,{\tau }_U)\mathbf{A}{({\mathit{\beta }}_0;{\tau }_L,{\tau }_U)}^{-1}$, where $\mathbf{A}({\mathit{\beta }}_0;{\tau }_L,{\tau }_U)$ is the limit of ${\mathbf{A}}_n({\mathit{\beta }}_0;{\tau }_L,{\tau }_U)$ and $\mathbf{V}({\mathit{\beta }}_0;{\tau }_L,{\tau }_U)$ is defined in the appendix.

In the appendix, we prove Theorems 1 and 2 by approximating $U_n(\mathit{\beta };{\tau }_L;{\tau }_U)$ with its smoothed analogy ${\tilde{U}}_n(\mathit{\beta };{\tau }_L,{\tau }_U)$ and hence showing its local asymptotic linearity in $\mathit{\beta }$. In the spirit of Ying (1993) and Lin and Ying (1995), the lack of monotonicity in rank-based estimating equations results in large sample properties in a local sense. Similar conclusion can also be seen in pointwise censored quantile regression models such as Peng and Huang (2008). Global consistency might be justified under particular special cases, see, for example, Jin et al. (2003) for a scenario where the accelerated failure time model with time-independent covariate and Gehan weight is considered.

4 Inference

4.1 Resampling Scheme

As shown in the proof of Theorem 2, the covariance matrix of the limiting distribution of $n^{1/2}(\widehat{\mathit{\beta }}-{\mathit{\beta }}_0)$ involves the unknown $\Lambda (t;{\mathit{\beta }}_0)$. To estimate the variances of the model parameters, we apply the resampling technique proposed in Jin, Ying, and Wei (2001) by perturbing the minimand of the estimating equations. Accordingly, we let $\{{\zeta }_1,\dots ,{\zeta }_n\}$ be iid nonnegative random variables with mean 1 and variance 1, such as exponential random variable with mean 1. We fix the data $\{{\tilde{T}}_i,{\Delta }_i,{\overline{\mathbf{Z}}}_i({\tilde{T}}_i),i=1,\dots ,n\}$ and generate $\{{\zeta }_1,\dots ,{\zeta }_n\}$, then consider the following perturbation of $U_n(\mathit{\beta };{\tau }_L,{\tau }_U)$, that is, $$U_n^{*}(\mathit{\beta };{\tau }_L,{\tau }_U)=\frac{1}{n}\sum _{i=1}^n{\zeta }_i\Psi (h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })$$ $$\times \left\{{\mathbf{Z}}_i({\tilde{T}}_i)-\frac{S^{(1)}(h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })}{S^{(0)}(h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })}\right\}\times$$ $$\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\times \{N_i(\widehat{q}({\tau }_U;\mathit{\beta });\mathit{\beta })-N_i(\widehat{q}({\tau }_L;\mathit{\beta });\mathit{\beta })\}.$$

Again we find the root of $U_n^{*}(\mathit{\beta };{\tau }_L,{\tau }_U)$, denoted by ${\mathit{\beta }}^{*}$, by minimizing its L₁ norm. Note that $\widehat{q}({\tau }_L;\mathit{\beta })$ and $\widehat{q}({\tau }_U;\mathit{\beta })$ are functions of $\mathit{\beta }$ which also depend on the perturbation as they are recomputed from the Kaplan-Meier estimator based on different values of the perturbed estimate ${\mathit{\beta }}^{*}$. While keeping the data fixed, we repeat the above procedure B times, that is, for each $r=1,\dots ,B$, generate a set of $\{{\zeta }_1,\dots ,{\zeta }_n\}$ and obtain a corresponding realization of ${\mathit{\beta }}^{*}$, denoted by ${\mathit{\beta }}_r^{*}$. Now, we state the following theorem, which will be proven in the appendix by invoking similar arguments as presented in Jin, Ying, and Wei (2001).

Theorem 3.

If Conditions 1–6 hold, then the conditional distribution of $n^{1/2}({\mathit{\beta }}^{*}-\widehat{\mathit{\beta }})$ given the observed data converges weakly to the same limiting distribution of $n^{1/2}(\widehat{\mathit{\beta }}-{\mathit{\beta }}_0)$.

As a consequence of Theorem 3, we may use $\{{\mathit{\beta }}_r^{*},r=1,\dots ,B\}$ to estimate the limiting distribution of the parameter estimates, that is, the variance of $\widehat{\mathit{\beta }}$ can be approximated by the sample variance of $\{{\mathit{\beta }}_r^{*},r=1,\dots ,B\}$, and the confidence interval for ${\mathit{\beta }}_0$ can be constructed using the empirical quantile of $\{{\mathit{\beta }}_r^{*},r=1,\dots ,B\}$ or by normal approximation.

4.2 Two-Stage Testing Procedure for Homogeneity

In this section, we present a procedure for testing the homogeneity across quantile regions. Based on the derivation of (2.4), the estimating equation remains valid as long as the homogeneity condition (2.1) holds over $[{\tau }_L,{\tau }_U]$. In such cases, the covariate effects represented by the regression coefficients are identical over the range of quantiles in which (2.1) holds. Hence, one may be interested in assessing the subregions of quantile over $[0,1]$, or more generally, any two quantile regions, for which the covariate effects are constant. This is equivalent to testing whether the regression coefficients over the subregions are identical or not.

For a specific interval $[{\tau }_L,{\tau }_U]\subset [0,1]$, suppose we want to test the constancy of ${\mathit{\beta }}_0$ within the range. A straightforward approach is to divide the interval into two halves, namely $[{\tau }_L,{\tau }_M]$ and $[{\tau }_M,{\tau }_U]$ where ${\tau }_M=({\tau }_L+{\tau }_U)/2$, and test whether or not the corresponding parameters in the two nonoverlapping regions are equal. More precisely, let ${\mathit{\beta }}_A$ denote the model parameter when the homogeneity condition (2.1) is assumed over a set A, and ${\widehat{\mathit{\beta }}}_A$ be the corresponding estimate using (2.4), it suffices to test if ${\mathit{\beta }}_{[{\tau }_L,{\tau }_M]}={\mathit{\beta }}_{[{\tau }_M,{\tau }_U]}$. One possible test is to reject the null hypothesis when the absolute deviation statistic $\mathcal{T}={\sum }_{k=1}^p|{\widehat{\mathit{\beta }}}_{[{\tau }_L,{\tau }_M]}^{(k)}-{\widehat{\mathit{\beta }}}_{[{\tau }_M,{\tau }_U]}^{(k)}|$ is large. As a result of Theorem 3, the rejection cutoff can be determined by the empirical quantile of the resampling-based ${\mathcal{T}}^{*}={\sum }_{k=1}^p|({\mathit{\beta }}_{[{\tau }_L,{\tau }_M]}^{*(k)}-{\mathit{\beta }}_{[{\tau }_M,{\tau }_U]}^{*(k)})-({\widehat{\mathit{\beta }}}_{[{\tau }_L,{\tau }_M]}^{(k)}-{\widehat{\mathit{\beta }}}_{[{\tau }_M,{\tau }_U]}^{(k)})|$. Theoretically, the test works for any arbitrary interval provided the sample size is large enough. However, because of the practical concern on the sample size of our data analysis, we propose a two stage testing procedure below that would identify the constant regions with reasonable statistical powers.

In particular, if one concerns about the lower quantiles $[0,0.25]$ of the failure time, she may question if the covariate effect is in fact homogeneous over $[0,1]$, and if not, whether it is same as that of $[0.25,0.5]$. Denote $A_1=[0,0.25],A_2=[0.25,0.5]$, $A_3=[0.5,0.75]$ and $A_4=[0.75,1]$, then the two questions translate to testing the two null hypotheses of $H_0^{(1)}:{\mathit{\beta }}_{A_1}={\mathit{\beta }}_{A_2}={\mathit{\beta }}_{A_3}={\mathit{\beta }}_{A_4}$ and $H_0^{(2)}:{\mathit{\beta }}_{A_1}={\mathit{\beta }}_{A_2}$, respectively. To this end, we first test $H_0^{(1)}$ using the test statistic ${\mathcal{T}}_1={\sum }_{k=1}^p\sum _{i\ne j}|{\widehat{\mathit{\beta }}}_{A_i}^{(k)}-{\widehat{\mathit{\beta }}}_{A_j}^{(k)}|$, where the rejection cutoff is taken to be the empirical quantile of ${\mathcal{T}}_1^{*}={\sum }_{k=1}^p\sum _{i\ne j}|({\mathit{\beta }}_{A_i}^{*(k)}-{\mathit{\beta }}_{A_j}^{*(k)})-({\widehat{\mathit{\beta }}}_{A_i}^{(k)}-{\widehat{\mathit{\beta }}}_{A_j}^{(k)})|$. Suppose $H_0^{(1)}$ is not rejected, we conclude that the covariate effects are identical across the subregions of quantile. On the other hand, if the first hypothesis $H_0^{(1)}$ is rejected, we proceed to test the second hypothesis of $H_0^{(2)}$. Similarly, we consider the test statistic ${\mathcal{T}}_2={\sum }_{k=1}^p|{\widehat{\mathit{\beta }}}_{A_1}^{(k)}-{\widehat{\mathit{\beta }}}_{A_2}^{(k)}|$ and use the empirical quantiles of ${\mathcal{T}}_2^{*}={\sum }_{k=1}^p|({\mathit{\beta }}_{A_1}^{*(k)}-{\mathit{\beta }}_{A_2}^{*(k)})-({\widehat{\mathit{\beta }}}_{A_1}^{(k)}-{\widehat{\mathit{\beta }}}_{A_2}^{(k)})|$ as cutoff. Rejection of $H_0^{(2)}$ means the covariate effects are different on A₁ and A₂.

The choices of A_i depend on the quantile range of interest. If the upper quantiles $[0.75,1]$ are concerned, one may set $A_1=[0.75,1],A_2=[0.5,0.75],A_3=[0.25,0.5]$ and $A_4=[0,0.25]$. Meanwhile, $A_1=[0.375,0.625]$, $A_2=[0.25,0.375]\cup [0.625,0.75],A_3=[0,0.25]$ and $A_4=[0.75,1]$ can be chosen if the median covariate effect is of interest. Despite the two-stage testing scheme, there is no need to adjust for the significance level for multiple hypotheses because of the nested structure of our hypotheses and procedure. Note that $H_0^{(1)}:{\mathit{\beta }}_{A_1}={\mathit{\beta }}_{A_2}={\mathit{\beta }}_{A_3}={\mathit{\beta }}_{A_4}$ is a subset of $H_0^{(2)}:{\mathit{\beta }}_{A_1}={\mathit{\beta }}_{A_2}$. Also, as we would only continue to test $H_0^{(2)}$ when we reject $H_0^{(1)}$, so rejecting $H_0^{(2)}$ is a subset of rejecting $H_0^{(1)}$. As a result, the probability of rejecting at least one hypothesis when both are true, also known as the family error rate, that is, $P(\text{Reject}{H}_0^{(1)}\cup \text{Reject}{H}_0^{(2)}\mid H_0^{(1)}\cap H_0^{(2)})$ reduces to $P(\text{Reject}{H}_0^{(1)}\mid H_0^{(1)})$. Recall that the family error rate is exactly the quantity that shall be controlled below the nominal significance level α to adjust for multiple tests (Simes 1986). Hence, we can set both $P(\text{Reject}{H}_0^{(1)}\mid H_0^{(1)})$ and $P(\text{Reject}{H}_0^{(2)}\mid H_0^{(2)})$ to be α, that is, the significance levels of both stages need not be modified, yet the Type I error of the procedure is controlled at α. Also, one may generalize the two-stage test by slicing $[0,1]$ into finer sub-regions and consider a multi-stage test so as to find out the region for which homogeneity holds more precisely, provided reasonable power can be achieved with the given sample size. Alternatively, the test can also serve as an a priori assessment of the suitable quantile interval upon the availability of a pilot study. Users may employ this test to validate their initial choice of quantile intervals for the full scaled study.

5 Simulations

In the following simulation examples, we consider a covariate process that changes in values exactly d – 1 times at the time nodes $\{W_{i1},\dots ,W_{i,d-1}\}$ where $0=W_{i0}<W_{i1}<\dots <W_{i,d-1}<W_{id}=\infty$. As a result, we can write $h_i(t;\mathit{\beta })$ as a piecewise linear function of t. Since $h_i(t;\mathit{\beta })$ is monotone in t, its inverse function $h_i^{-1}(t;\mathit{\beta })$ is also piecewise linear in t. Note that by the equivariance property of the quantile function, model (2.1) can be written as $Q_{T_i}(u\mid {\overline{\mathbf{Z}}}_i)=h_i^{-1}(q(u);{\mathit{\beta }}_0)$ for $u\in [{\tau }_L,{\tau }_U]$, which can be exploited to generate the desired lifetime with a designated conditional quantile function by the inverse transformation method.

In particular, let ${\mathbf{Z}}_i(t)={(Z_{i1}(t),Z_{i2}(t))}^{\top }$ be a two-dimensional covariate process with d = 3 subintervals. Specifically, let $$Z_{i1}(t)=\sum _{s=1}^d{X}_{i1s}I(W_{i,s-1}\le t<W_{is}),$$ $$Z_{i2}(t)=\sum _{s=1}^d{X}_{i2s}I(W_{i,s-1}\le t<W_{is}),$$where $X_{i1s},X_{i2s}\sim \text{Gamma}(2,1/2)$ independently for each $i=1,\dots ,n$ and $s=1,\dots ,d$. The true regression coefficient is ${\mathit{\beta }}_0={(1,1)}^{\top }$ and the common conditional quantile function $\hspace{0.17em}\log \hspace{0.17em}q(\tau )=Q_{\eta }(\tau )$ where $\eta \sim N(0,1)$. For $i=1,\dots ,n$ and $s=1,\dots ,d-1$, the time nodes are recovered recursively through $W_{is}=W_{i,s-1}+S_{is}$, where S_is (the sojourn times between the time nodes) are independently distributed as U(0, 1). Recall that the covariate processes $\{{\mathbf{Z}}_i(s),0\le s\le {\tilde{T}}_i\}$ can only be observed up to the survival times, so we generate the sojourn time in a way that the simulated failure times $\{T_i,i=1,\dots ,n\}$ spread through the subintervals derived from $\{W_{i0},\dots ,W_{id}\}$. Although we consider two time-dependent covariates in the simulation studies, any time-independent covariates can be easily incorporated by treating them as constant processes. One can also notice that the resulting conditional quantile function of the failure time remains monotone increasing for any realization of the covariate process.

Example 1:

Homogeneous Error

Given the above setup, the first numerical example demonstrates that the proposed method resembles the classic accelerated failure time model when the homogeneous conditional quantile assumption (2.1) is satisfied for $[{\tau }_L,{\tau }_U]=[0,1]$. To this end, we generate the censoring times $\{C_i,i=1,\dots ,n\}$ from iid exponential distribution with mean $\hspace{0.17em}\exp \hspace{0.17em}\{c+100Z_{i1}(0)+100Z_{i2}(0)\}$ so that the censoring times depend on the covariate processes, where we set c = – 150 or c = – 110 to achieve a censoring rate of 40% or 20%, respectively. With ${\mathit{\beta }}_0={(1,1)}^{\top }$ and $\hspace{0.17em}\log \hspace{0.17em}q(\tau )=Q_{\eta }(\tau )$, we have (2.1) holds for $[{\tau }_L,{\tau }_U]=[0,1]$, and thus for a subinterval $[{\tau }_L,{\tau }_U]=[0,0.5]$. We use (2.4) with $[{\tau }_L,{\tau }_U]=[0,0.5]$ to obtain our parameter estimates $\widehat{\mathit{\beta }}$. We generate 500 Monte Carlo datasets with n = 200 or n = 400 sample points each, and adopt the aforementioned resampling procedure with B = 200 and ζ_i exponentially distributed with mean 1 to reconcile the asymptotic distribution of the parameter estimates. The numerical optimization is carried out via the differential evolution algorithm in MATLAB. The choice of differential evolution parameters follows the rules of thumb suggested by Storn (1996). Further discussion on the parameter selection can also be found in Price, Storn, and Lampinen (2006). We adopt existing quantile regression techniques with time-independent covariates, such as Wang and Wang’s (2009) approach, to obtain preliminary initial estimates. The average value over time of the covariate process observed in the study is considered as a time-independent proxy to be used in their methods. Our numerical experience suggests that the monotonicity of Wang and Wang’s (2009) estimating equation with respect to the parameter of interest serves as a reasonable guideline for initial value selection.

Remark 1.

It is noteworthy that differential evolution is not the only viable algorithm for the optimization problem concerned. Indeed, as suggested by an anonymous referee, we may update $\widehat{\mathit{\beta }}$ and $\widehat{q}$ in (2.4) in an alternate manner, that is, update $\widehat{q}$ given the current $\widehat{\mathit{\beta }}$ and solve for $\widehat{\mathit{\beta }}$ by fixing $\widehat{q}$ iteratively, until convergence. Simulation results obtained by this iterative method is presented in the appendix for reference. We observe that this method achieves similar numerical results as compared to the differential evolution algorithm, with slightly larger empirical biases under some circumstances. We trust that the iterative method can be a valuable and effective alternative especially for scenarios where each update of $\widehat{\mathit{\beta }}$ involves a convex and smooth optimization problem. Unfortunately, this simplification may not bring in significant improvements for the problem discussed in this article or nonconvex optimization problems in general.

Table 1 reports the simulation results for our proposal in terms of empirical bias (Bias), empirical standard deviation (SD), average standard error (SE) based on the resampling method, and coverage rate (CP) of the 95% confidence intervals constructed by normal approximation. The results obtained by our proposal using the different weight functions are quite similar, so the estimation is not sensitive with respect to the specific choice of weight function. The empirical performance is compared against Lin and Ying’s (1995) and Gorfine, Goldberg, and Ritov’s (2017) proposals. Under the current setting where the accelerated failure time model holds, the proposed method produces competitive results in empirical biases and standard deviations when compared to Lin and Ying (1995). Moreover, the standard errors are closed to the empirical standard deviations and the coverage probabilities are close to the nominal level of 95%, validating the aforementioned resampling method. Our proposal loses some efficiency against Lin and Ying’s (1995) method as reflected in slightly larger empirical standard deviations. This is reasonable because we only consider (2.4) over a subinterval $[{\tau }_L,{\tau }_U]=[0,0.5]$ in which (2.1) holds to obtain $\widehat{\mathit{\beta }}$, while they consider the whole region $[{\tau }_L,{\tau }_U]=[0,1]$. On the other hand, despite Gorfine, Goldberg, and Ritov’s (2017) method only requires (2.1) to hold for a single value of τ, their estimates demonstrate larger empirical biases and standard deviations, especially when the censoring rate is high. One reason is that their methodology only works for cases with independent censoring, which obviously does not hold in this example as the censoring distribution depends on the covariate processes.

Table 1 Simulation results with homogeneous error.

		n = 200				n = 400
Censoring rate		20%		40%		20%		40%
		${\widehat{\beta }}^{(1)}$	${\widehat{\beta }}^{(2)}$	${\widehat{\beta }}^{(1)}$	${\widehat{\beta }}^{(2)}$	${\widehat{\beta }}^{(1)}$	${\widehat{\beta }}^{(2)}$	${\widehat{\beta }}^{(1)}$	${\widehat{\beta }}^{(2)}$
Proposed (unweighted)	Bias	−0.009	−0.011	−0.008	−0.011	0.000	−0.007	−0.001	−0.008
	SD	0.149	0.150	0.162	0.171	0.105	0.105	0.121	0.118
	SE	0.151	0.152	0.168	0.168	0.105	0.105	0.119	0.117
	CP	0.946	0.942	0.954	0.938	0.940	0.936	0.952	0.940
Proposed (Gehan weighted)	Bias	−0.009	−0.009	−0.010	−0.013	0.000	−0.007	0.000	−0.006
	SD	0.146	0.147	0.160	0.167	0.103	0.104	0.119	0.116
	SE	0.146	0.148	0.163	0.165	0.103	0.102	0.115	0.115
	CP	0.938	0.952	0.954	0.946	0.946	0.938	0.944	0.938
Lin and Ying (1995)	Bias	−0.006	−0.009	−0.002	−0.011	−0.004	−0.005	−0.004	−0.004
	SD	0.132	0.127	0.149	0.146	0.092	0.093	0.108	0.104
Gorfine, Goldberg, and Ritov (2017) ($\tau =0.25$)	Bias	−0.767	−0.763	−0.849	−0.817	−0.722	−0.717	−0.755	−0.732
	SD	0.959	1.078	0.701	0.743	0.564	0.507	0.408	0.685
Gorfine, Goldberg, and Ritov (2017) ($\tau =0.50$)	Bias	−0.768	−0.807	−0.910	−0.783	−0.700	−0.743	−0.739	−0.817
	SD	1.082	1.901	1.452	1.540	1.154	0.421	0.764	1.484
Gorfine, Goldberg, and Ritov (2017) ($\tau =0.75$)	Bias	−0.748	−0.706	−0.766	−0.529	−0.789	−0.804	−0.794	−0.760
	SD	2.029	2.251	1.869	6.726	0.391	0.453	0.364	0.613
Augmented—Gorfine, Goldberg, and Ritov (2017) ($\tau =0.25$)	Bias	−0.102	−0.070	−0.139	−0.133	−0.132	−0.136	−0.058	−0.113
	SD	0.425	0.414	0.404	0.411	0.394	0.377	0.473	0.411
Augmented—Gorfine, Goldberg, and Ritov (2017) ($\tau =0.50$)	Bias	−0.055	−0.057	−0.085	−0.065	−0.063	−0.062	−0.065	−0.095
	SD	0.322	0.319	0.363	0.361	0.357	0.332	0.423	0.393
Augmented—Gorfine, Goldberg, and Ritov (2017) ($\tau =0.75$)	Bias	−0.041	−0.041	−0.052	−0.048	−0.038	−0.043	−0.059	−0.058
	SD	0.301	0.308	0.322	0.327	0.303	0.308	0.311	0.351
Composite—Gorfine, Goldberg, and Ritov (2017)	Bias	−0.839	−0.847	−0.898	−0.890	−0.770	−0.751	−0.774	−0.780
	SD	0.983	1.157	0.721	0.770	0.604	0.530	0.426	0.722

We also employ the augmentation-based estimator discussed in Gorfine, Goldberg, and Ritov (2017) for a more thorough comparison, which we denote as Augmented- Gorfine, Goldberg, and Ritov (2017) in Tables 1–3. The augmentation component exploits the information from the censored observations to yield a doubly-robust estimator. To extend its scope so that it can handle conditional independent censoring, we incorporate a kernel-based local Kaplan-Meier estimator (see, Wang and Wang 2009) for estimating the conditional censoring distribution given the covariates into the augmentation-based estimator. Equipped with a more refined estimate for the censoring distribution, the augmented estimator manages to outperform its vanilla counterpart in terms of empirical bias. However, as observed in Gorfine, Goldberg, and Ritov (2017), the augmented version may not necessarily improve the empirical standard deviations. Our proposed method still maintains a reasonable edge, particularly on the empirical standard deviations. It is noteworthy that the local Kaplan-Meier estimator is kernel-based and requires the use of instrumental variables for the time-dependent covariates. However, the performance of the arbitrary choice of the instrumental variables is fairly sensitive to the censoring mechanism. In practice, as the underlying censoring mechanism can hardly be determined, this casts concerns in the actual implementation. Another drawback of the augmentation-based estimator as described by the authors is the long computing time compared to its vanilla version and our proposal. In our numerical experience with a E5-2650v4 CPU and 8GB memory, the augmented method requires about 21 min to run for a single quantile, against 1 sec of its vanilla counterpart and 20 sec of the proposed estimator.

Table 2 Simulation results with homogeneity assumed for $\tau \le \mathrm{0.25.}$

		n = 200				n = 400
Censoring rate		20%		40%		20%		40%
		${\widehat{\beta }}^{(1)}$	${\widehat{\beta }}^{(2)}$	${\widehat{\beta }}^{(1)}$	${\widehat{\beta }}^{(2)}$	${\widehat{\beta }}^{(1)}$	${\widehat{\beta }}^{(2)}$	${\widehat{\beta }}^{(1)}$	${\widehat{\beta }}^{(2)}$
Proposed (unweighted)	Bias	−0.017	−0.014	−0.023	−0.023	−0.001	−0.007	−0.001	−0.005
	SD	0.201	0.203	0.238	0.230	0.134	0.136	0.163	0.156
	SE	0.198	0.201	0.233	0.233	0.137	0.135	0.154	0.155
	CP	0.944	0.952	0.946	0.950	0.948	0.934	0.926	0.954
Proposed (Gehan weighted)	Bias	−0.016	−0.013	−0.017	−0.016	0.000	−0.005	0.001	−0.000
	SD	0.201	0.202	0.230	0.222	0.135	0.134	0.155	0.150
	SE	0.196	0.198	0.227	0.229	0.134	0.133	0.152	0.150
	CP	0.950	0.950	0.954	0.946	0.944	0.932	0.930	0.946
Lin and Ying (1995)	Bias	0.809	0.804	0.810	0.796	0.822	0.820	0.819	0.820
	SD	0.183	0.181	0.213	0.220	0.129	0.130	0.155	0.145
Gorfine, Goldberg, and Ritov (2017) ($\tau =0.25$)	Bias	−0.609	−0.650	−0.685	−0.753	−0.495	−0.495	−0.596	−0.592
	SD	0.626	0.633	0.516	0.492	0.427	0.431	0.409	0.418
Gorfine, Goldberg, and Ritov (2017) ($\tau =0.50$)	Bias	−0.753	−0.802	−0.845	−0.863	−0.735	−0.741	−0.790	−0.786
	SD	0.614	0.548	0.533	0.571	0.513	0.514	0.531	0.482
Gorfine, Goldberg, and Ritov (2017) ($\tau =0.75$)	Bias	−1.701	−1.733	−1.750	−1.782	−1.445	−1.493	−1.421	−1.415
	SD	0.769	0.722	0.670	0.677	0.772	0.812	0.838	0.828
Augmented—Gorfine, Goldberg, and Ritov (2017) ($\tau =0.25$)	Bias	0.029	−0.016	−0.021	0.018	0.159	0.188	0.233	0.300
	SD	0.639	0.681	0.688	0.785	0.813	0.747	0.939	0.958
Augmented—Gorfine, Goldberg, and Ritov (2017) ($\tau =0.50$)	Bias	0.225	0.154	0.141	0.141	0.258	0.263	0.119	0.174
	SD	0.930	0.944	0.995	0.915	1.076	0.990	0.948	1.016
Augmented—Gorfine, Goldberg, and Ritov (2017) ($\tau =0.75$)	Bias	0.005	0.034	−0.022	−0.041	0.054	0.033	0.019	0.000
	SD	0.372	0.370	0.437	0.431	0.347	0.337	0.411	0.406
Composite—Gorfine, Goldberg, and Ritov (2017)	Bias	−0.688	−0.727	−0.758	−0.811	−0.551	−0.551	−0.627	−0.659
	SD	0.740	0.733	0.535	0.575	0.484	0.559	0.476	0.488

Table 3 Simulation results with homogeneity assumed for $0.25\le \tau \le \mathrm{0.75.}$

		n = 200				n = 400
Censoring rate		20%		40%		20%		40%
		${\widehat{\beta }}^{(1)}$	${\widehat{\beta }}^{(2)}$	${\widehat{\beta }}^{(1)}$	${\widehat{\beta }}^{(2)}$	${\widehat{\beta }}^{(1)}$	${\widehat{\beta }}^{(2)}$	${\widehat{\beta }}^{(1)}$	${\widehat{\beta }}^{(2)}$
Proposed (unweighted)	Bias	0.079	0.079	0.081	0.053	0.047	0.061	0.058	0.068
	SD	0.341	0.348	0.433	0.424	0.230	0.239	0.295	0.276
	SE	0.314	0.315	0.375	0.379	0.210	0.215	0.254	0.256
	CP	0.930	0.940	0.922	0.946	0.958	0.942	0.930	0.948
Proposed (Gehan weighted)	Bias	0.028	0.015	0.006	−0.019	0.016	0.016	0.005	0.015
	SD	0.227	0.218	0.318	0.315	0.123	0.134	0.208	0.157
	SE	0.199	0.199	0.261	0.262	0.118	0.121	0.152	0.151
	CP	0.960	0.966	0.958	0.958	0.950	0.956	0.946	0.968
Lin and Ying (1995)	Bias	0.228	0.203	0.241	0.206	0.192	0.205	0.193	0.209
	SD	0.363	0.358	0.404	0.404	0.237	0.238	0.271	0.255
Gorfine, Goldberg, and Ritov (2017) ($\tau =0.25$)	Bias	−1.318	−1.350	−1.594	−1.676	−1.039	−1.048	−1.419	−1.366
	SD	0.811	0.781	0.855	0.860	0.716	0.752	0.644	0.662
Gorfine, Goldberg, and Ritov (2017) ($\tau =0.50$)	Bias	−0.827	−0.805	−0.924	−0.922	−0.499	−0.493	−0.600	−0.585
	SD	0.721	0.726	0.899	0.989	0.578	0.605	0.589	0.543
Gorfine, Goldberg, and Ritov (2017) ($\tau =0.75$)	Bias	−0.804	−0.740	−0.844	−0.870	−0.582	−0.592	−0.683	−0.663
	SD	0.606	0.656	0.618	0.588	0.516	0.507	0.458	0.473
Augmented—Gorfine, Goldberg, and Ritov (2017) ($\tau =0.25$)	Bias	−0.690	−0.695	−1.139	−1.222	−0.742	−0.723	−1.156	−1.157
	SD	0.637	0.655	0.638	0.588	0.568	0.601	0.585	0.501
Augmented—Gorfine, Goldberg, and Ritov (2017) ($\tau =0.50$)	Bias	−0.113	−0.076	−0.179	−0.175	0.074	−0.052	−0.064	−0.086
	SD	0.515	0.521	0.572	0.598	0.637	0.565	0.732	0.669
Augmented—Gorfine, Goldberg, and Ritov (2017) ($\tau =0.75$)	Bias	0.052	0.135	−0.008	0.128	0.119	−0.010	−0.110	0.007
	SD	0.889	0.914	0.907	0.968	0.941	0.876	0.840	1.000
Composite—Gorfine, Goldberg, and Ritov (2017)	Bias	−0.837	−0.826	−0.898	−0.956	−0.511	−0.495	−0.589	−0.589
	SD	0.697	0.676	0.725	0.782	0.614	0.607	0.630	0.601

The composite quantile regression introduced in Zou and Yuan (2008) also shares the pooling effect that our methodology provides, although the two ideas have different modeling motivations. We include their possible extension to censored data for comparison, which we refer to as Composite- Gorfine, Goldberg, and Ritov (2017). The modification exploits and composites the most relevant existing quantile regression model for time-dependent covariates, that is, Gorfine, Goldberg, and Ritov’s (2017) methodology, over a grid of quantiles. We use a grid size of 0.1 over the interested range of quantiles $[0,0.5]$ for illustration here. While the composite feature could possibly improves the numerical performance of the ordinary quantile regression, it involves a higher dimensional optimization procedure, hence, the numerical problem turns out to be more challenging here. It is shown that the Gorfine, Goldberg, and Ritov’s (2017) method shares similar empirical bias with Gorfine, Goldberg, and Ritov (2017) while its standard deviations are slightly inferior, which are both not comparable to the proposed method. The main reason is that the estimates being composited requires an independent censoring assumption which does not hold in our examples.

Example 2:

Heterogeneity in Upper Tail Quantiles

The second example is more illustrative to show that our proposed method generalizes the accelerated failure time model proposed by Lin and Ying (1995), in the sense that it relaxes homogeneous conditional quantile assumption (2.1) for all quantile levels. With the same covariate processes described above, we now generate the failure times based on $Q_{T_i}(\tau \mid {\overline{\mathbf{Z}}}_i)=h_i^{-1}(q(\tau );{\mathit{\beta }}_0+I(\tau >0.5))$. Here, ${\mathit{\beta }}_0+I(\tau >0.5)$ reduces to ${\mathit{\beta }}_0$ for $\tau \le 0.5$, which implies that model (2.1) holds for $[{\tau }_L,{\tau }_U]=[0,0.5]$. But the expression becomes ${\mathit{\beta }}_0+1$, not the same constant ${\mathit{\beta }}_0$, for $[{\tau }_L,{\tau }_U]=[0.5,1]$, therefore, model (2.1) is not valid over the entire $[{\tau }_L,{\tau }_U]=[0,1]$. Owing to the heterogeneity across the two subregions of quantiles, the accelerated failure time model assumed by Lin and Ying (1995) does not hold here. Nevertheless, suppose we are only interested in the covariate effect on the lower quantiles $[0,0.25]$, then our proposal of using (2.4) with $[{\tau }_L,{\tau }_U]=[0,0.25]$ to compute $\widehat{\mathit{\beta }}$ will still be a justified methodology. Under the same simulation setting as in the first example, Table 2 reports the simulation results in the same format as Table 1 does. It is clear that the proposed estimator produces reasonable results as in the homogeneous setting described in the first example, because our model (2.1) holds for $[{\tau }_L,{\tau }_U]=[0,0.25]$ here. Meanwhile, it is reasonable that Lin and Ying’s (1995) method produces substantially larger empirical bias and standard deviation than it does in the first example, because of the violation of their model assumption. The performance of vanilla, augmented, and composite Gorfine, Goldberg, and Ritov’s (2017) methods are worse than our proposal as in the first example, which is due to the same reason of the violation of the independent censoring assumption.

Example 3:

Heterogeneity in Both Lower and Upper Tail Quantiles

The third example considers failure times with quantile function $$Q_{T_i}(\tau \mid {\overline{\mathbf{Z}}}_i)=h_i^{-1}(q(\tau );{\mathit{\beta }}_0-Q_{\eta }(\tau )I(\tau <0.25)$$ $$+ Q_{\eta }(\tau )I(\tau >0.75)).$$

The second argument of the last display reduces to ${\mathit{\beta }}_0$ for $\tau \in [0.25,0.75]$. Therefore, model (2.1) holds for $[{\tau }_L,{\tau }_U]=[0.25,0.75]$ whereas such equivalence does not hold beyond that. Suppose we want to examine the covariate effect on the failure time around the median, we may consider a subinterval of quantiles around the median, say $[{\tau }_L,{\tau }_U]=0.5\pm 0.25=[0.25,0.75]$ for the estimating equation (2.4)(2.4) $$\begin{array}{c}{U}_n(\mathit{\beta };{\tau }_L,{\tau }_U)=\frac{1}{n}\sum _{i=1}^n\Psi (h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })\\ \times \left\{{\mathbf{Z}}_i({\tilde{T}}_i)-\frac{S^{(1)}(h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })}{S^{(0)}(h_i({\tilde{T}}_i;\mathit{\beta });\mathit{\beta })}\right\}\\ \hspace{1em}\times \{N_i(\widehat{q}({\tau }_U;\mathit{\beta });\mathit{\beta })-N_i(\widehat{q}({\tau }_L;\mathit{\beta });\mathit{\beta })\}\mathrm{}.\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\end{array}$$(2.4) . As shown in Table 3, the empirical biases and standard deviations of our method can mostly demonstrate the advantage as compared to the competitors. Same as the previous scenario, Lin and Ying’s (1995) method, Gorfine, Goldberg, and Ritov’s (2017) vanilla proposal, augmented method, and composite modification yield significantly larger biases because of violation of the homogeneous conditional quantile assumption over all quantiles and the independent censoring assumption, respectively.

With a significance level of 5%, Table 4 reports the empirical rejection rate of the two-stage test described in Section 4 for the previous three simulation examples with the unweighted estimating equation. We assume $A_1=[0,0.25],A_2=[0.25,0.5],A_3=[0.5,0.75]$ and $A_4=[0.75,1]$ for the first two examples; $A_1=[0.375,0.625],A_2=[0.25,0.375]\cup [0.625,0.75],A_3=[0,0.25]$ and $A_4=[0.75,1]$ for the third example. It is shown that the test has well controlled empirical Type I errors with reasonable powers under the alternative hypothesis.

Table 4 Empirical rejection rate of the test for homogeneity.

	n = 200		n = 400
Censoring rate	20%	40%	20%	40%
Situations where the null is true
Example 1, ${\mathcal{T}}_1$	0.040	0.030	0.028	0.030
Example 1, ${\mathcal{T}}_2$	0.034	0.028	0.052	0.050
Example 2, ${\mathcal{T}}_2$	0.048	0.052	0.050	0.070
Example 3, ${\mathcal{T}}_2$	0.008	0.020	0.018	0.026
Situations where the null is false
Example 2, ${\mathcal{T}}_1$	0.904	0.806	0.988	0.976
Example 3, ${\mathcal{T}}_1$	0.632	0.492	0.778	0.644

Example 4:

Independent Censoring

Despite the proposed method assumes a local homogeneity around the interested quantile level compared to the quantile regression method described by Gorfine, Goldberg, and Ritov (2017), it effectively exploits more information than the latter to allow a more numerically stable estimation, provided such an assumption holds. Here, we present a simulation example to demonstrate that the stability in the proposed method offers a more reliable inference across a range of quantiles even under the independent censoring setup. We consider covariate processes as described in the previous examples, while the true regression coefficients are now ${\mathit{\beta }}_0={(0.8,0.8)}^{\top }$. The failure times are generated based on the conditional quantile function of $Q_{T_i}(\tau \mid {\overline{\mathbf{Z}}}_i)=h_i^{-1}(q(\tau );{\mathit{\beta }}_0)$, while the censoring times are simulated from iid exponential distribution with mean (3.6) to yield a censoring rate of 20%. The failure times and the censoring times are then marginally independent in which case the condition assumed in Gorfine, Goldberg, and Ritov (2017) holds. Their method is implemented for every quantile level from 0.25 to 0.75 with an interval of 0.05, while the proposed method is adopted with $[{\tau }_L,{\tau }_U]=[0.25,0.75]$ for comparison.

Table 5 shows that when the independent censoring assumption holds, Gorfine, Goldberg, and Ritov’s (2017) proposal provides reasonable estimates over the quantile levels considered. Nevertheless, their empirical confidence intervals for both parameter estimates are volatile and fluctuate across neighboring quantiles. This phenomenon is not desirable as confusing conclusions may be drawn on whether or not a covariate is significant over the inter-quantile range concerned. As a hypothetical example, if the user is interested in the failure time distribution around the median, looking at $\tau =0.5$ would conclude that the first covariate is insignificant while the second is significant. But if one performs the analysis by checking the neighboring quantile $\tau =0.45$ or $\tau =0.55$, a different conclusion will be drawn, which may be confusing. Here, we imagine analyzing the designated quantile level separately, so there is no multiple comparison involved. On the other hand, provided that the homogeneity assumption holds over $[{\tau }_L,{\tau }_U]$, our proposal pools information over this range to produce a more numerically stable estimate, so that the conclusion can be more consistently drawn. We also study the performance of the Composite- Gorfine, Goldberg, and Ritov (2017) variant discussed above for reference. This composite method produces estimates that are averages of the vanilla Gorfine, Goldberg, and Ritov’s (2017) estimates. As compared to the proposed method, they have larger empirical biases and standard deviations.

Table 5 Simulation results with homogeneous error and independent censoring.

	${\widehat{\beta }}^{(1)}$		${\widehat{\beta }}^{(2)}$
	Estimate	95% CI	Estimate	95% CI
Proposed ($[{\tau }_L,{\tau }_U]=[0.25,0.75]$)	0.777	$( 0.444, 1.110)$	0.776	$( 0.450, 1.103)$
Gorfine, Goldberg, and Ritov (2017)
$\tau =0.25$	0.771	$(-0.007, 1.549)$	0.804	$(-0.042, 1.649)$
$\tau =0.30$	0.749	$(-0.011, 1.508)$	0.814	$(-0.025, 1.652)$
$\tau =0.35$	0.788	$( 0.021, 1.555)$	0.772	$(-0.031, 1.576)$
$\tau =0.40$	0.740	$(-0.016, 1.496)$	0.747	$(-0.014, 1.509)$
$\tau =0.45$	0.736	$( 0.049, 1.422)$	0.734	$( 0.070, 1.399)$
$\tau =0.50$	0.731	$(-0.015, 1.476)$	0.754	$( 0.030, 1.478)$
$\tau =0.55$	0.725	$(-0.015, 1.464)$	0.710	$(-0.005, 1.426)$
$\tau =0.60$	0.694	$( 0.011, 1.377)$	0.699	$(-0.079, 1.476)$
$\tau =0.65$	0.672	$( 0.005, 1.339)$	0.676	$( 0.050, 1.302)$
$\tau =0.70$	0.677	$(-0.004, 1.358)$	0.647	$(-0.010, 1.305)$
$\tau =0.75$	0.658	$(-0.003, 1.318)$	0.649	$(-0.083, 1.381)$
Composite - Gorfine, Goldberg, and Ritov (2017)	0.726	$( 0.008, 1.444)$	0.750	$(-0.008, 1.509)$

6 Data Analysis

In this section, we demonstrate the proposed method with the Stanford heart transplant data. The study is described in details in Clark et al. (1971) and Crowley and Hu (1977); see also Aitkin, Laird, and Francis (1983) for a brief literature review on the statistical analysis of the data. Lin and Ying (1995) analyzed the effect of heart transplantation on the lifetime using the accelerated failure time model with time-dependent covariates. Their results suggest that younger patients with lower mismatch scores are more likely to benefit from transplantation, while it is less beneficial for older patients with higher mismatch scores to receive transplant. Here, in particular, we are interested in examining the effect of heart transplantation in a neighborhood of the median lifetime, and those in the lower quantile and the upper quantile.

Same as Lin and Ying (1995), we consider a three-dimensional time-dependent covariate process in our analysis, namely transplant status, age at transplant, and mismatch score, where the last two covariates are assumed to be zero before transplant. Mathematically, let W_i denote the waiting time to transplant of the ith patient, that is, the duration from the date of acceptance into the program to the date of transplant. We take W_i to be $\infty$ if the ith patient never received transplantation. If we denote ${\mathbf{Z}}_i(t)={\{Z_{i1}(t),Z_{i2}(t),Z_{i3}(t)\}}^{\top }$ to be the three-dimensional covariate process of the ith patient, then $$\begin{array}{cc}{Z}_{i1}(t)\hfill & =I(t\ge W_i),\hfill \\ Z_{i2}(t)\hfill & =(\text{age at transplant}-35)I(t\ge W_i),\hfill \\ Z_{i3}(t)\hfill & =(\text{mismatch score}-0.5)I(t\ge W_i).\hfill \end{array}$$

Here, we centralize $Z_{i2}(t)$ and $Z_{i3}(t)$ so that the regression coefficient associated with $Z_{i1}(t)$ can be interpreted as the effect of transplantation for a patient who is 35 years old at the time of transplant and has a mismatch score of 0.5.

Suppose we are interested in the covariate effects on the bottom quantiles, neighborhood of median, and the top quantiles. We may set $[{\tau }_L,{\tau }_U]$ to be $[0,0.25],[0.375,0.625]$ and $[0.5,1]$, respectively. Here, $[0.5,1]$ is considered instead of $[0.75,1]$ to avoid any possible identifiability issues on the extreme quantiles. Table 6 reports the parameter estimates and their corresponding 95% confidence intervals. The 95% confidence intervals are constructed by normal approximation under the aforementioned resampling scheme with B = 500. Estimation results based on the accelerated failure time model, the Cox model, and the quantile model (at $\tau =0.25,0.5,0.75$), each under the framework of time-dependent covariates, are also provided as in Gorfine, Goldberg, and Ritov (2017). With a significance level of 5%, both the accelerated failure time model and the Cox model conclude that the transplant status and the age at transplant are significant, whereas the mismatch score is not. Based on Gorfine, Goldberg, and Ritov’s (2017) approach, the quantile regression models at 25th, 50th, and 75th percentiles reveal that the covariate effects are not significant with the exception of transplant status at 25th and 50th percentiles, which might be due to the higher variance associated with their estimator as shown in the previous simulation examples. While the composite variant of Gorfine, Goldberg, and Ritov (2017) generally produces smaller standard errors than the vanilla version, its estimates for the bottom quantiles are relatively unstable. Possible causes can be the higher dimensional optimization problem involved in the composite method and the fact that Gorfine, Goldberg, and Ritov’s (2017) procedure is intrinsically volatile across extreme tail quantiles in nature. On the other hand, the parameter estimates of the proposed method largely share consistent signs with those given by Lin and Ying (1995) and the Cox model. The estimated coefficients for transplant status are larger on the middle and top quantiles than that on bottom quantile range, suggesting that a transplantation can be more beneficial to those patients who live longer than their counterparts with similar covariate background. Meanwhile, the age at transplant has a negative influence on the failure time over all three regions of quantiles. For the mismatch score, the estimated coefficients flip from positive on the bottom quantile range to negative on the higher quantile ranges. It can be interpreted as a transplant mismatch would be more unfavorable to the patients surviving longer after covariate adjustment, which again infers that a transplantation is more influential to that cohort. Such observations could be reasonable in nature, because it would take time for the transplantation to become effective. Comparing the estimated coefficients for mismatch score across the three regions, the difference can be an indication that the mismatch score having an insignificant overall effect over all quantile levels as estimated by Lin and Ying (1995).

Table 6 Parameter estimates of the Stanford heart transplant data.

		Transplant status	Age at transplant–35	Mismatch score–0.5
Proposed
$[{\tau }_L,{\tau }_U]=[0,0.25]$	Estimate	0.379	–0.035	$0.322$
	95% CI	$(-0.585, 1.343)$	$(-0.124, 0.054)$	$(-0.534, 1.178)$
$[{\tau }_L,{\tau }_U]=[0.375,0.625]$	Estimate	4.134	–0.083	–2.352
	95% CI	$( 3.285, 4.983)$	$(-0.144,-0.021)$	$(-2.896,-1.809)$
$[{\tau }_L,{\tau }_U]=[0.5,1]$	Estimate	2.067	–0.015	–1.322
	95% CI	$( 1.908, 2.227)$	$(-0.028,-0.001)$	$(-1.527,-1.116)$
Gorfine, Goldberg, and Ritov (2017)
$\tau =0.25$	Estimate	2.553	$0.083$	–0.708
	95% CI	$( 0.885, 4.220)$	$(-0.310, 0.475)$	$(-2.654, 1.237)$
$\tau =0.5$	Estimate	2.452	–0.062	–0.513
	95% CI	$( 0.637, 4.267)$	$(-0.238, 0.114)$	$(-1.371, 0.346)$
$\tau =0.75$	Estimate	2.013	–0.009	–1.196
	95% CI	$(-0.399, 4.425)$	$(-0.455, 0.437)$	$(-2.466, 0.739)$
Composite—Gorfine, Goldberg, and Ritov (2017)
${\tau }_j\in [0,0.25]$	Estimate	$1.181$	–0.243	$0.701$
	95% CI	$(-3.769, 6.132)$	$(-1.224, 0.738)$	$(-3.515, 4.916)$
${\tau }_j\in [0.375,0.625]$	Estimate	$2.016$	0.032	–1.220
	95% CI	$( 1.794, 2.238)$	$(-0.044, 0.108)$	$(-1.572,-0.869)$
${\tau }_j\in [0.5,1]$	Estimate	$4.299$	$0.208$	–0.438
	95% CI	$( 3.647, 4.951)$	$(-0.097, 0.513)$	$(-0.956, 0.081)$
Lin and Ying (1995)	Estimate	1.986	–0.096	–0.930
	p-value	0.028	$0.003$	$0.155$
Cox Model	Estimate	1.031	–0.055	–0.445
	p-value	0.033	$0.015$	$0.112$

We also implement the two-stage testing procedure described in Section 4 for assessing the quantile regions of homogeneous covariate effects under a 5% level of significance. With $[{\tau }_L,{\tau }_U]=[0,0.25]$, the empirical p-value of the first stage is 0.01, so rejecting $H_0^{(1)}:{\mathit{\beta }}_{A_1}={\mathit{\beta }}_{A_2}={\mathit{\beta }}_{A_3}={\mathit{\beta }}_{A_4}$ suggests that the error is not homogeneous. Furthermore, the second stage produces an empirical p-value of 0.01, thus, rejecting $H_0^{(2)}:{\mathit{\beta }}_{A_1}={\mathit{\beta }}_{A_2}$ implies the covariate effects over $[0,0.25]$ and $[0.25,0.5]$ are different. For $[{\tau }_L,{\tau }_U]=[0.375,0.625]$, the empirical p-values of the two stages are 0 and 0, respectively, that is, both $H_0^{(1)}$ and $H_0^{(2)}$ are rejected, hence, the homogeneous quantile region concluded would be $[0.375,0.625]$. Finally, for $[{\tau }_L,{\tau }_U]=[0.75,1]$, the empirical p-values are 0.01 and 0.224, so $H_0^{(2)}$ is not rejected in the second stage, that is, the covariate effects over $[0.75,1]$ and $[0.5,0.75]$ are not distinguishable. If further numerical study is to be carried out on a full scaled dataset, this result strengthens the idea of using $[0.5,1]$ to access the covariate effects on the top quantiles.

7 Discussion

In view of this new perspective in survival analysis, further studies may also be of interest for generalizing the proposed methodology to accommodate different variations of failure time data in practice. Particularly, modification can be made without disrupting the framework of rank-based estimation procedure so as to handle length-biased data and case-cohort studies. On the other hand, following the exploration of quantile localized covariate effects, one may abandon the assumption of universal covariate effects but assume a hierarchical structure of the regression parameters instead. Being viewed as a hybrid between our setup and the censored quantile regression model, the current proposal can be generalized when the regression parameters are presumed to be with certain structure, especially provided prior knowledge on the covariate effects. Investigation of these new directions shall be handled in separate papers.

Supplementary Materials

The supplementary materials contain additional details about the asymptotic results and the simulation studies. Section A to C of the appendix present respectively the proofs of Theorems 1–3 discussed in Sections 3 and 4. Section D summarizes the numerical results for the alternative iterative algorithm discussed in Remark 1 of Section 5.

Acknowledgments

The authors would like to thank the editor, the associate editor and the anonymous referees for their constructive comments and valuable suggestions with which the manuscript have been substantially improved.

Additional information

Funding

Sit’s work was partially supported by Hong Kong Research Grant Council RGC-14301618, RGC-14301920, and RGC-14307221.