Nonparametric quantile regression for twice censored data

Nonparametric quantile regression for twice censored data

Stanislav Volgushev, Holger Dette,
Ruhr-Universität Bochum
University of Illinois at Urbana-Champaign.
Supported by the Sonderforschungsbereich “Statistical modeling of nonlinear dynamic processes” (SFB 823) of the Deutsche Forschungsgemeinschaft.

We consider the problem of nonparametric quantile regression for twice censored data. Two new estimates are presented, which are constructed by applying concepts of monotone rearrangements to estimates of the conditional distribution function. The proposed methods avoid the problem of crossing quantile curves. Weak uniform consistency and weak convergence is established for both estimates and their finite sample properties are investigated by means of a simulation study. As a by-product, we obtain a new result regarding the weak convergence of the Beran estimator for right censored data on the maximal possible domain, which is of its own interest.

AMS Subject Classification: 62G08, 62N02, 62E20

Keywords and Phrases: quantile regression, crossing quantile curves, censored data, monotone rearrangements, survival analysis, Beran estimator

1 Introduction

Quantile regression offers great flexibility in assessing covariate effects on event times. The method was introduced by Koenker and Bassett (1978) as a supplement to least squares methods focussing on the estimation of the conditional mean function and since this seminal work it has found numerous applications in different fields [see Koenker (2005)]. Recently Koenker and Geling (2001) have proposed quantile regression techniques as an alternative to the classical Cox model for analyzing survival times. These authors argued that quantile regression methods offer an interesting alternative, in particular if there is heteroscedasticity in the data or inhomogeneity in the population, which is a common phenomenon in survival analysis [see Portnoy (2003)]. Unfortunately the “classical” quantile regression techniques cannot be directly extended to survival analysis, because for the estimation of a quantile one has to estimate the censoring distribution for each observation. As a consequence rather stringent assumptions are required in censored regression settings. Early work by Powell (1984, 1986), requires that the censoring times are always observed. Moreover, even under this rather restrictive and – in many cases – not realistic assumption the objective function is not convex, which results in some computational problems [see for example Fitzenberger (1997)]. Even worse, recent research indicates that using the information contained in the observed censored data actually reduces the estimation accuracy [see Koenker (2008)].
Because in most survival settings the information regarding the censoring times is incomplete several authors have tried to address this problem by making restrictive assumptions on the censoring mechanism. For example, Ying et al. (1995) assumed that the responses and censoring times are independent, which is stronger than the usual assumption of conditional independence. Yang (1999) proposed a method for median regression under the assumption of i.i.d. errors, which is computationally difficult to evaluate and cannot be directly generalized to the heteroscedastic case. Recently, Portnoy (2003) suggested a recursively re-weighted quantile regression estimate under the assumption that the censoring times and responses are independent conditionally on the predictor. This estimate adopts the principle of self consistency for the Kaplan-Meier statistic [see Efron (1967)] and can be considered as a direct generalization of this classical estimate in survival analysis. Peng and Huang (2008) pointed out that the large sample properties of this recursively defined estimate are still not completely understood and proposed an alternative approach, which is based on martingale estimating equations. In particular, they proved consistency and asymptotic normality of their estimate.

While all of the cited literature considers the classical linear quantile regression model with right censoring, less results are available for quantile regression in a nonparametric context. Some results on nonparametric quantile regression when no censoring is present can be found in Chaudhuri (1991) and Yu and Jones (1997, 1998). Chernozhukov et al. (2006) and Dette and Volgushev (2008) pointed out that many of the commonly proposed parametric or nonparametric estimates lead to possibly crossing quantile curves and modified some of these estimates to avoid this problem. Results regarding the estimation of the conditional distribution function from right censored data can be found in Dabrowska (1987, 1989) or Li and Doss (1995). The estimation of conditional quantile functions in the same setting is briefly stressed in Dabrowska (1987) and further elaborated in Dabrowska (1992a), while El Ghouch and Van Keilegom (2008) proposed a quantile regression procedure for right censored and dependent data. On the other hand, the problem of nonparametric quantile regression for censored data where the observations can be censored from either left or right does not seem to have been considered in the literature.

This gap can partially be explained by the difficulties arising in the estimation of the conditional distribution function with two-sided censored data. The problem of estimating the (unconditional) distribution function for data that may be censored from above and below has been considered by several authors. For an early reference see Turnbull (1974). More recent references are Chang and Yang (1987); Chang (1990); Gu and Zhang (1993) and Patilea and Rolin (2006). On the other hand- to their best knowledge- the authors are not aware of literature on nonparametric conditional quantile regression, or estimation of a conditional distribution function, for left and right censored data when the censoring is not always observed and only the conditional independence of censoring and lifetime variables is assumed.

In the present paper we consider the problem of nonparametric quantile regression for twice censored data. We consider a censoring mechanism introduced by Patilea and Rolin (2006) and propose an estimate of the conditional distribution function in several steps. On the basis of this estimate and the preliminary statistics which are used for its definition, we construct two quantile regression estimates using the concept of simultaneous inversion and isotonization [see Dette et al. (2005)] and monotone rearrangements [see Dette et al. (2006), Chernozhukov et al. (2006) or Anevski and Fougères (2007) among others]. In Section 2 we introduce the model and the two estimates, while Section 3 contains our main results. In particular, we prove uniform consistency and weak convergence of the estimates of the conditional distribution function and its quantile function. As a by-product we obtain a new result on the weak convergence of the Beran estimator on the maximal possible interval, which is of independent interest. In Section 4 we illustrate the finite sample properties of the proposed estimates by means of a simulation study. Finally, all proofs and technical details are deferred to an Appendix.

2 Model and estimates

We consider independent identically distributed random vectors , , where are the variables of interest, and are left and right censoring variables, respectively, and the -valued random variables denote the covariates. We assume that the distributions of the random variables and depend on and denote by the conditional distribution function of given . The conditional distribution functions and are defined analogously.

Additionally, we assume that the random variables are almost surely nonnegative and independent conditionally on the covariate . Our aim is to estimate the conditional quantile function . However, due to the censoring, we can only observe the triples where and the indicator variables are defined by

Remark 2.1

An unconditional version of this censoring mechanism was introduced by Patilea and Rolin (2006). Examples of situations where this kinds of data occur can for example be found in chapter of Meeker and Escobar (1998). This model also is closely related to the double censoring model, see Turnbull (1974) for the case without covariates. In that setting, the assumption of independence between the random variables is replaced by the assumption that is independent of the pair and additionally . Note that none of the two assumptions is strictly more or less restrictive then the other. Rather the two models describe different situations. Moreover, since are never observed simultaneously, it is not possible to decide which of the models is most approriate. Instead, an understanding of the underlying data generation process is crucial to identify the right model. A more detailed comparison of the two models can be found in Patilea and Rolin (2001) and Patilea and Rolin (2006) for the case without covariates.

Roughly speaking, the construction of an estimate for the conditional quantile function of can be accomplished in three steps. First, we define the variables and consider the model , which is a classical right censoring model. In this model we estimate the conditional distribution of . In a second step, we use this information to reconstruct the conditional distribution of [see Section 2.1]. Finally, the concept of simultaneous isotonization and inversion [see Dette et al. (2005)] and the monotone rearrangements, which was recently introduced by Dette et al. (2006) in the context of monotone estimation of a regression function, are used to obtain two estimates of the conditional quantile function [see Section 2.2].

2.1 Estimation of the conditional distribution function

To be more precise, let denote the conditional distribution of . We introduce the notation and obtain the decomposition for the conditional distribution of . The sub-distribution functions can be represented as follows


Note that the conditional (sub-)distribution functions and can easily be estimated from the observed data by


where the quantities denote local weights depending on the covariates , which will be specified below. We will use the representations (2.2) - (2.4) to obtain an expression for in terms of the functions and then replace the distribution functions by their empirical counterparts , respectively. We begin with the reconstruction of . First note that


is the predictable reverse hazard measure corresponding to and hence we can reconstruct using the product-limit representation


[see e.g. Patilea and Rolin (2006)]. Now having a representation for the conditional distribution function we can define in a second step

which yields an expression for the predictable hazard measure of . Finally, can be reconstructed by using the product-limit representation


[see e.g. Gill and Johansen (1990)]. Note that formula (2.9) yields an explicit representation of the conditional distribution function in terms of the quantities , which can be estimated from the data [see equation (2.5)]. The estimate of the conditional distribution function is now defined as follows. First, we use the representation (2.7) to obtain an estimate of , that is




Second, after observing (2.1) and (2.9), we define




In Section 3 we will analyse the asymptotic properties of these estimates, while in the following Section 2.2 these estimates are used to construct nonparametric and noncrossing quantile curve estimates.

Remark 2.2

Throughout this paper, we will adopt the convention . This means that if, for example, and , the contribution of

in (2.13) will be interpreted as zero.

2.2 Non-crossing quantile estimates by monotone rearrangements

In practice, nonparametric estimators of a conditional distribution function are not necessarily increasing for finite sample sizes [see e.g. Yu, Jones (1998)]. Although this problem often vanishes asymptotically, it still is of great practical relevance, because in a concrete application it is not completely obvious how to invert a non-increasing function. Trying to naively invert such estimators may lead to the well-known problem of quantile crossing [see Koenker (2005) or Yu and Jones (1998)] which poses some difficulties in the interpretation of the results. In this paper we will discuss the following two possibilities to deal with this problem

  1. Use a procedure developed by Dette and Volgushev (2008) which is based on a simultaneous isotononization and inversion of a nonincreasing distribution function. As a by-product this method yields non-crossing quantile estimates. To be precise, we consider the operator


    where denotes the set of bounded, measurable functions on the set and denotes a bounded interval. Note that for a strictly increasing function this operator yields the right continuous inverse of , that is [here and in what follows, will denote the generalized inverse, i.e. ]. On the other hand, is always isotone, even in the case where does not have this property. Consequently, if is a not necessarily isotone estimate of an isotone function , the function could be regarded as an isotone estimate of the function . Therefore, the first idea to construct an estimate of the conditional quantile function consists in the application of the operator to the estimate defined in (2.12), i.e.


    However, note that formally the mapping operates on functions defined on bounded intervals. More care is necessary if the operator has to be applied to a function with an unbounded support. A detailed discussion and a solution of this problem can be found in Dette and Volgushev (2008). In the present paper we use different approach which is a slightly modified version of the ideas from Anevski and Fougères (2007). To be precise note that estimators of the conditional distribution function [in particular those of the form (2.5), which will be used later] often are constant outside of the compact interval . Now the structure of the estimator implies that will also be constant outside of . We thus propose to consider the modified operator defined as


    Consequently the first estimator of the conditional quantile function is given by

  2. Use the concept of increasing rearrangements [see Dette et al. (2006) and Chernozhukov et al. (2006) for details] to construct an increasing estimate of the conditional distribution function, which is then inverted in a second step. More precisely, we define the operator


    where is introduced in (2.14). Note that for a strictly increasing right continuous function this operator reproduces , i.e. . On the other hand, if is not isotone, is an isotone function and the operator preserves the -norm, i.e.

    Moreover, the operator also defines a contraction, i.e.

    [see Hardy et al. (1988) or Lorentz (1953)]. This means if is a not necessarily isotone estimate of the isotone function , then the isotonized estimate is a better approximation of the isotone function than the original estimate with respect to any -norm [note that because is assumed to be isotone]. For a general discussion of monotone rearrangements and the operators (2.14) and (2.18) we refer to Bennett and Sharpley (1988), while some statistical applications can be found in Dette et al. (2006) and Chernozhukov et al. (2006).
    The idea is now to use rearranged estimators of and in the representations (2.6)-(2.9). For this purpose we need to modify the operator so that it can be applied to functions of unbounded support. We propose to proceed as follows

    • Define the operator indexed by the compact interval as

    • Truncate the estimator for values outside of the interval , i.e.

      [note that in general estimators of the form (2.5) do not necessarily have values in the interval since the weights might be negative]

    • Use the statistic as estimator for .

    • Observe that the estimator is by construction an increasing step function which can only jump in the points , i.e. it admits the representation


      with weights . Based on this statistic, we define estimators of the subdistribution functions as follows


      In particular, such a definition ensures that .

    So far we have obtained increasing estimators of the quantities and . The next step in our construction is to plug these estimates in representation (2.6) to obtain:


    which defines an increasing function with jumps of size less or equal to one. This implies that is also increasing. For the rest of the construction, observe the following Lemma which will be proved at the end of this section.

    Lemma 2.3

    Assume that for . Then the function


    is nonnegative, increasing and has jumps of size less or equal to one.

    This in turn yields the estimate


    In the final step we now simply invert the resulting estimate of the conditional distribution function since it is increasing by construction. We denote this estimator of the conditional quantile function by


In the next section, we will discuss asymptotic properties of the two proposed estimates and of the conditional quantile curve.

Remark 2.4

In the classical right censoring case, there is no uniformly good way to define the Kaplan-Meier estimator beyond the largest uncensored observation [see e.g. Fleming and Harrington (1991), page 105]. Typical approaches include setting it to unity, to the value at the largest uncensored observation, or to consider it unobservable within certain bounds [for more details, see the discussion in Fleming and Harrington (1991), page 105 and Anderson et al. (1993), page 260]. When censoring is light, the first of the above mentioned approaches seems to yield the best results [see Anderson et al. (1993), page 260].
When the data can be censored from either left or right, the situation becomes even more complicated since now we also have to find a reasonable definition below the smallest uncensored observation. From definitions (2.6)-(2.9) it is easy to see that equals zero below the smallest uncensored observation with non-vanishing weight and is constant at the largest uncensored observation and above. In practice, the latter implies that the estimators and are not defined as soon as or , respectively. A simple ad-hoc solution to this problem is to define the estimator or as beyond the last observation with non-vanishing weight or to locally increase the bandwidth. A detailed investigation of this problem is postponed to future research.

We conclude this section with the proof of Lemma 2.3.
Proof of Lemma 2.3 In order to see that is increasing, we note that

Thus and the nonnegativity of is established. In order to prove the inequality we assume without loss of generality that . Observe that as soon as we have for

where the equalities and follow from . An analogous result for follows by simple algebra. Hence we have established that for we have , and all the other cases need not be considered since we adopted the convention ’0/0=0’. Thus the proof is complete.

3 Main results

The results stated in this section describe the asymptotic properties of the proposed estimators. In particular, we investigate weak convergence of the processes , etc. where the predictor is fixed. Our main results deal with the weak uniform consistency and the weak convergence of the process and the corresponding quantile processes obtained in Section 2. In order to derive the process convergence, we will assume that it holds for the initial estimates and give sufficient conditions for this property in Lemma 3.3. In a next step we apply the delta method [see Gill (1989)] to the map defined in and the product-limit maps defined in and . Note that the product limit maps are Hadamard differentiable on the set of cadlag functions with total variation bounded by a constant [see Lemma A.1 on page 42 in Patilea and Rolin (2001)], and hence the process convergence of and will directly entail the weak convergence results for and , respectively. However, the Hadamard differentiability of the map only holds on domains where , and hence more work is necessary to obtain the corresponding weak convergence results on the interval if , where


This situation occurs for example if , which is quite natural in the context considered in this paper because is the right censoring variable.

For the sake of a clear representation and for later reference, we present all required technical conditions for the asymptotic results at the beginning of this section. We assume that the estimators of the conditional subdistribution functions are of the form (2.5) with weights depending on the covariates but not on or . The first set of conditions concerns the weights that are used in the representation (2.5). Throughout this paper, denote by the maximum norm on .

  1. [label=(W0)]

  2. With probability tending to one, the weights in (2.5) can be written in the form

    where the real-valued functions have the following properties:

    1. [label=(0), ref=1(0)]

    2. There exist constants such that for all and all we have either or

    3. If for some constant , then and for for some sequence such that . Without loss of generality, we will assume that throughout this paper.

    4. for some positive function .

    5. .

    Here [and throughout this paper] denotes a smoothing parameter converging to 0 with increasing sample size.

  3. We assume that the weak convergence

    holds in , where the limit denotes a centered Gaussian process which has a version with a.s. continuous sample paths and a covariance structure of the form

    for some function . Here and throughout this paper weak convergence is understood as convergence with respect to the sigma algebra generated by the closed balls in the supremum norm [see Pollard (1984)].

  4. The estimators and are weakly uniformly consistent on the interval

Remark 3.1

It will be shown in Lemma 3.3 below that, under suitable assumptions on the smoothing parameter , important examples for weights satisfying conditions 1-3 are given by the Nadaraya-Watson weights


or (in one dimension) by the local linear weights

where , and the kernel satisfies the following condition.

  1. The kernel in (3.2) and (3.1) is a symmetric density of bounded total variation with compact support, say , which satisfies for all with for some constants .

For the distributions of the random variables we assume that for some with

  1. [label=(D0)]

  2. The conditional distribution function fulfills

  3. For we have

  4. The conditional distribution functions have densities,
    say , with respect to the Lebesque measure

  5. The functions are twice continuously differentiable with respect to the second component in some neighborhood of and for we have

  6. The distribution function of the covariates is twice continuously differentiable in . Moreover, has a uniformly continuous density such with .

  7. There exists a constant such that for all where is a set with the property for some and all .

  8. uniformly in as

  9. For we have

Remark 3.2

From the definition of and we immediately see that under condition 1 we have where we use the notation . In particular, this implies that under either of the assumptions 4 or 11 the equality holds.

Finally, we make some assumptions for the smoothing parameter

  1. [label=(B0)]

  2. and .

  3. and .

Some important practical examples for weights satisfying conditions 1 - 3 include Nadaraya-Watson and local linear weights. This is the assertion of the next Lemma.

Lemma 3.3
  1. Conditions 1a and 1b are fulfilled for the Nadaraya-Watson weights with a Kernel satisfying condition (K1). If the density is continuous at the point , condition 1c also holds. Finally, if the function is continuously differentiable in a neighborhood of for every with uniformly (in ) bounded first derivative and (B1) is fulfilled, condition 1d holds.
    If additionally to these assumptions and the density of the covariates is continuously differentiable at with bounded derivative, condition 1 also holds for the local linear and rearranged local linear weights and defined in (3.1) and (2.20), (2.21) respectively, provided that the corresponding kernel fulfills condition (K1) .

  2. If under assumptions 7, 8 and 1 the density is twice continuously differentiable with uniformly bounded derivative, condition 2 holds for the Nadaraya-Watson ( arbitrary), local linear () or rearranged local linear () weights based on a positive, symmetric kernel with compact support.

  3. If under assumptions 2, 2, 3 the density is twice continuously differentiable with uniformly bounded derivative, condition 3 holds for the Nadaraya-Watson weights based on a positive, symmetric kernel with compact support ( arbitrary). If additionally and the density of the covariates is continuously differentiable at with bounded derivative, condition 3 also holds for local linear or rearranged local linear weights.

The proof of this Lemma is standard, a sketch can be found in the Appendix.

Note that the assumption 1 does not allow to choose , which would be the MSE-optimal rate for Nadaraya-Watson or local linear weights and functions with two continuous derivatives with respect to the predictor. This assumption has been made for the sake of a transparent presentation and implies that the bias of the estimates is negligible compared to the stochastic part. Such an approach is standard in nonparametric estimation for censored data, see Dabrowska (1987) or Li and Doss (1995). In principle, most results of the present paper can be extended to bandwidths if a corresponding bias term is subtracted.

Another useful property of estimators constructed from weights satisfying condition 1 is that they are increasing with probability tending to one.

Lemma 3.4

Under condition 1a we have

The Lemma follows from the relation

and the fact that under assumption (W1) the probability of the event on the right hand side converges to one. We will use Lemma 3.4 for the analysis of the asymptotic properties of the conditional quantile estimators in Section 3.2. One noteworthy consequence of the Lemma is the fact that

which follows because the mappings and the right continuous inversion mapping coincide on the set of nondecreasing functions. In particular, this indicates that, from an asymptotic point of view, it does not matter which of the estimators is used. The difference between both estimators will only be visible in finite samples - see Section 4. In fact, it can only occur if one of the estimators is decreasing at some point.

3.1 Weak convergence of the estimate of the conditional distribution

We are now ready to describe the asymptotic properties of the estimates defined in Section 2. Our first result deals with the weak uniform consistency of the estimate under some rather weak conditions. In particular, it does neither require the existence of densities of the conditional distribution functions [see 3] nor integrability conditions like 4.

Theorem 3.5

If conditions 1, 3, 11, 1a-1b and 3 are satisfied, then the following statements are correct.

  1. The estimate defined in (2.12) is weakly uniformly consistent on the interval for any such that .

  2. If additionally