Semiparametric Copula Quantile Regression for Complete or Censored Data

When facing multivariate covariates, general semiparametric regression techniques come at hand to propose flexible models that are unexposed to the curse of dimensionality. In this work a semiparametric copula-based estimator for conditional quantiles is investigated for complete or right-censored data. In spirit, the methodology is extending the recent work of Noh et al. (2013) and Noh et al. (2015), as the main idea consists in appropriately defining the quantile regression in terms of a multivariate copula and marginal distributions. Prior estimation of the latter and simple plug-in lead to an easily implementable estimator expressed, for both contexts with or without censoring, as a weighted quantile of the observed response variable. In addition, and contrary to the initial suggestion in the literature, a semiparametric estimation scheme for the multivariate copula density is studied, motivated by the possible shortcomings of a purely parametric approach and driven by the regression context. The resulting quantile regression estimator has the valuable property of being automatically monotonic across quantile levels, and asymptotic normality for both complete and censored data is obtained under classical regularity conditions. Finally, numerical examples as well as a real data application are used to illustrate the validity and finite sample performance of the proposed procedure.


Introduction
Quantile regression is a prevailing method when it comes to investigating the possible relationships between a d-dimensional covariate X and a response variable T . Since the seminal work of Koenker and Basset (1978), quantile regression has received notable interest in the literature on theoretical and applied statistics as a very attractive alternative to the classical mean regression model based on quadratic loss. As the latter only captures the central tendancy of the data, there are many cases and nice examples where mean regression is uninformative with respect to studying the conditional upper or lower quantiles. For an interesting application, see for example Elsner et al. (2008). A comprehensive review of quantile regression as a robust (to outliers) and flexible (to error distribution) method can be found in Koenker (2005).
A wide literature on the estimation of the quantile regression function is devoted to the case where the response variable T is completely observed. This is not necessarily the case in survival analysis, where right censoring of T may arise, i.e. instead of fully observing the variable of interest, one only observes the minimum of it and a censoring variable. For instance, in clinical studies, censoring may occur because of the withdrawal of patients from the study, the end of the follow-up period, etc. In this context, quantile regression becomes attractive as an alternative to popular regression techniques like the Cox proportional hazards model or the accelerated failure time model, as is argued in Koenker and Bilias (2001), Koenker and Geling (2001) and Portnoy (2003). Additional appealing properties of the method include the fact that it allows for modelling heterogeneity of the variance and it does not necessarily impose a proportional effect of the covariates on the hazard over the duration time as opposed to the popular Cox model.
As is the case for the uncensored situation, existing literature on censored quantile regression includes fully parametric, semiparametric and nonparametric methodologies. When several covariates are to be taken into account, fully parametric methodologies are known to be highly sensitive to model misspecification and may lack the flexibility needed for an adequate modelling. On the other hand, in spite of their great flexibility, fully nonparametric methods such as local linear smoothing proposed by El Ghouch and Van Keilegom (2009) are typically affected by the curse of dimensionality. In light of these restrictions, semiparametric estimation procedures such as a single-index model suggested by Bücher et al. (2014) come at hand when the dimension of the covariate is high. It is the object of this paper to extend the rather sparse literature on flexible multidimensional methodologies in the context of censored quantile regression.
Censored quantile regression was first introduced by Powell (1986) for linear models and fixed censoring, that is, presuming that the censoring times are the same for all observations. For random censoring, Ying et al. (1995) proposed a semiparametric linear median regression model, assuming that the survival time and the censoring variable are unconditionally independent. Despite this important contribution, the proposed procedure involves solving non-monotone discontinuous equations, hereby introducing practical and computational difficulties. Furthermore, the unconditional independence assumption may sound restrictive as conditional independence, given the covariates, seems more natural in certain applications.
Based on conditional independance between the survival time and the censoring variable, Portnoy (2003) developed a novel estimating procedure based on the idea of redistribution-of-mass introduced by Efron (1967). A major shortcoming of it is the need for a global linear assumption, that is, in order to estimate the τ -th conditional quantile, one needs to assume the linearity of all the conditional functionals at lower quantiles. Extending Portnoy's work, Wang and Wang (2009) relaxed this constraint by only assuming linearity at one pre-specified quantile level of interest. Recently, Leng and Tong (2013) proposed an alternative to Wang and Wang's methodology by extending the work of Ying et al. in order to relax the constraining unconditional independance assumption and provide an efficient algorithm for estimation. Still, a linear approach may be too restrictive for real data applications.
Finally, an interesting alternative approach to the above-mentioned literature was proposed by Bücher et al. (2014), where a single-index model for the conditional quantile function is stud-ied under the assumption of independence between the covariates and the censoring variable. The single-index structure assumes that the objective function depends linearly on the covariates through an unknown link function, making the proposed model (i.e. under the aforementioned assumption) insensitive to the curse of dimensionality since the nonparametric part is of dimension one.
In this paper, we aim to extend the literature on multivariate quantile regression estimation in the possible presence of censored data by providing a rich, flexible and robust alternative based on the copula function that defines the dependence structure between the variables of interest. By taking advantage of copula modelling, we intend to provide a new class of estimators that would allow practitioners to analyse, in a flexible way, multidimensional survival data. Actually, our methodology is, in essence, an extension of the recent work of Noh et al. (2013) and Noh et al. (2015), as the central idea is to express the conditional quantile function in terms of an appropriate copula density and marginal distributions. In their original paper, Noh et al. (2013) suggested subsequently to leave the marginal distributions unspecified while assuming a parametric model for the copula. Overall, their suggested approach results in a semiparametric regression estimator that is not exposed to the curse of dimensionality. However, in order to avoid possible shortcomings highlighted by Dette et al. (2014) that are induced by the misspecification of the parametric copula, we propose in this work, both for complete and censored data, an alternative semiparametric estimation strategy for the copula itself. Our resulting methodology is flexible for multidimensional data with or without censoring, easy to implement and does not require any iterative procedure in opposition to existing semiparametric alternatives.
The rest of this paper is organised as follows. Developing the copula-based estimation procedure for the quantile regression is the topic of Section 2. The asymptotic properties of the proposed estimator are obtained in Section 3 and the finite sample performance is illustred by means of Monte Carlo simulations in Section 4, where both the semiparametric copula estimation strategy and the overall performance of our estimator are investigated. Section 5 provides a brief application to real data. Lastly, the proofs of our asymptotic properties are deferred to the Appendix.

Copula-based estimator for complete data
Let X = (X 1 , . . . , X d ) T be a covariate vector of dimension d ≥ 1 and T be a (time-to-event) response variable with marginal continuous cumulative distribution functions (c.d.f.) F 1 , . . . , F d and F T , respectively. Throughout this paper, we denote by f j and f T the density of X j , j = 1, . . . , d, and T , respectively. From the pioneering work of Sklar (1959), for a given x = (x 1 , . . . , x d ) T , the c.d.f. of (T, X) evaluated at (t, x) can be expressed as From Sklar's theorem, it is clear that the copula C T X disjoints the marginal behaviours of T and X from their dependence structure, hence allowing a great modelling flexibility. For a book length treatment of copulas, see Nelsen (2006) and Joe (2014).
The object of interest of this paper, the τ -th conditional quantile function of the dependent variable T given X = x, denoted by m τ (x), is defined for any τ ∈ (0, 1) as m τ (x) = inf{t : where ρ τ (u) = u(τ − 1(u < 0)) is the so-called "check" function, and 1(·) is the indicator function.
To motivate our approach, let us suppose, for the moment, that there is no censoring and that we observe an i.i.d. sample (T i , X i ), i = 1, . . . , n, from (T, X). In this context, following the definition of a copula function, Noh et al. (2015) noted that the conditional quantile function of T given X = x may be expressed as . . ∂u d is the copula density corresponding to C T X . Consequently, any given estimators p c T X , p F T and p F j of c T X , F T and F j , j = 1, . . . , d, respectively, automatically yield an estimator of m τ (x) given by As indicated earlier, Noh et al. suggest to estimate the marginals nonparametrically and to consider a parametrization of the copula density, that is, assume that the latter belongs to a certain parametric family of copula densities C = {c(u 0 , u; θ), θ ∈ Θ ⊂ R p }.
In this paper however, we propose a novel semiparametric strategy for the estimation of the copula density, motivated by the issues related to the possible misspecification of the parametric approach. To highlight this shortcoming and illustrate how one may circumvent it, we consider the simplistic example reported by Dette et al. (2014) with a single covariate, where (T i , X i ), i = 1, . . . , n, are i.i.d. random variables with T i = (X i − 0.5) 2 + σ i , X i ∼ U [0, 1], σ = 0.025 and i , i = 1, . . . , n, are i.i.d. standard normal random variables. In this situation, where the true quantile regression function is non-monotonic in the covariate, it is found that most of the common parametric copula families still yield a monotone estimation of the regression function, thereby providing a rather poor fit of the latter. This is illustrated in Figure 1a, where the estimation is carried out for τ = 0.5, n = 500 and using three common parametric copulas.
As the roots of the above-mentioned limitation are not intrinsic to a copula-based approach, but rather to be attributed to the limited set of parametric copula families existing in the literature, a natural alternative, for low dimensional covariates only, would be to consider a fully nonparametric estimation of the copula density itself. The resulting, and adequate, quantile regression estimation is depicted in Figure 1b.
Recalling that we intend to handle multivariate covariates in this paper, we will not adopt a purely nonparametric approach, but rather prefer a copula estimation strategy that provides sufficient flexibility to the multidimensional estimator while dodging dimension related constraints. More specifically, we note that any multivariate copula density can be decomposed into two parts as follows: where F j|T , j = 1, . . . , d, denotes the conditional c.d.f. of X j given T . The first part of the decomposition contains the product of d bivariate copula densities related to the dependence of T with every covariate, whereas the second part captures the conditional dependence of X given T = t.
In the general regression context, part (2.4) may then be interpreted as the dependence of actual interest since it focuses on the relationship between the response variable with every covariate. On the contrary, part (2.5) may be viewed in such framework as a 'noisy' dependence, or, more precisely, a correction parameter for possible (conditional) dependence among covariates. Consequently, a natural reasoning suggests to provide as much flexibility as possible to the modelling of part (2.4), while keeping the estimation of part (2.5) uncomplicated. We therefore advocate to estimate nonparametrically the d bivariate copulas of interest and, subsequently, exploit standard parametric techniques for the second part of the multivariate copula density.
Concerning the nonparametric estimation of bivariate copula densities, several methodologies have been proposed in the literature. To estimate for instance c T X 1 , a tempting, yet naive, approach would be to opt for standard multivariate kernel techniques. That is, given a bivariate sample (U 0i , U 1i ), i = 1, . . . , n, from (F T (T ), F 1 (X 1 )), one would estimate the copula density c T X 1 at where H is a symmetric positive-definite bandwidth matrix and K : R 2 → R is a bivariate kernel function. This technique is called naive in this context since it completely ignores the fact that the density of interest is only supported on the unit square. This boundedness property is at the origin of distinct schemes designed to correct for well-known bias issues at the boundaries, see for example the mirror reflection method (Gijbels and Mielniczuk (1990)) or the boundary kernel method (Chen and Huang (2007)). An appealing alternative was furthermore proposed by Charpentier et al. (2006) and Geenens et al. (2014), where the main idea is to appropriately project the initial data on an unbounded support with the purpose of then estimating the obtained bivariate transformed density by means of standard techniques (standard kernel (Charpentier et al.) or polynomial locallikelihood (Geenens et al.)). Using the invariance property of copulas to increasing transformations of their margins, the estimation of the copula density is then obtained by back-transformation on the unit square. That is, using the example of a probit transformation, one may estimate the copula density c T X 1 at (u 0 , u 1 ) ∈ [0, 1] 2 by where φ and Φ stand for the standard normal density and c.d.f., respectively, and where p f 01 is a bivariate density estimator of the projected data (Φ −1 (U 0i ), Φ −1 (U 1i )), i = 1, . . . , n. This transformation technique, coupled with polynomial local-likelihood estimation for f 01 , in order to allow for possible unbounded copula density estimates, is shown to outperform its competitors in most scenarios in a detailed simulation study in Geenens et al.. An exhaustive comparison of existing methodologies for bivariate estimation may be found in Nagler (2014). Furthermore, fully nonparametric multidimensional copulas are studied in Hobaek Haff and Segers (2012) and Nagler and Czado (2015).
Regarding the second part (2.5) of the (d + 1)-dimensional copula, standard high-dimensional parametric techniques involve, among others, nested Archimedean copulas (see e.g. Hofert and Pham (2013) and Joe (2014)), factor copulas (see e.g. Oh and Patton (2012)) and, arguably the most popular, vine copulas (see e.g. Czado (2010), Joe (2014) and references therein). Note, however, that for the estimation of (2.5) in practice, one is first advised to adopt the so-called simplifying assumption which stipulates that the conditioning on T = t is fully captured by the conditional marginals. In other words, in (2.5), the conditional copula itself is not affected by the conditioning on T . This assumption turns out to be the cornerstone of vine copula models as it keeps them tractable for inference and model selection. For more details about this and its implications, see Hobaek Haff et al. (2010) and Stöber et al. (2013).
Conclusively, in this article we propose to adopt the following detailed procedure for the modelling and estimation of the multivariate copula density: (1) Based on original observations (T i , X 1i , . . . , X di ), i = 1, . . . , n, construct 'pseudo-observations', needed for the estimation of (2.4), using rescaled versions of empirical distributions: where the factor 1/(n + 1), commonly adopted in the copula literature, aims at keeping the constructed observations in the interior of [0, 1].
(3) Compute the pseudo-observations needed for the estimation of (2.5) as: This relationship is at the origin of the sequential nature of the vine copula estimation scheme (see e.g. Czado (2010)).

Copula-based estimator for censored data
In the presence of censoring, the estimation equation (2.3) becomes inappropriate as we do not fully observe the response variables T i . Instead, we only observe a sequence of i.i.d. triplets (Y i , X i , ∆ i ), i = 1, . . . , n, from (Y, ∆, X), where Y = min(T, C), ∆ = 1(T ≤ C) and C denotes the censoring variable, assumed to be independent of T given X. In order to take censoring into account in the estimation procedure, the first step is to note that, for any measurable function ϕ : R → R, where G C (c|x) = P(C ≤ c|X = x) denotes the conditional distribution of C given X = x. This, along with (2.1), suggests a natural way to handle censoring for quantile regression by replacing the function ϕ with the check function. At this stage, a naive way of trying to take profit of copula modelling would be to consider introducing copulas in the obtained conditional expectation. Note, however, that the conditional expectation in (2.6) is in fact the joint conditional expectation of (Y, ∆) given X = x. Adopting an analogous reasoning as the one presented by Noh et al. for the uncensored case at this point would therefore result in the insertion of the joint copula of (Y, ∆, X), hence exposing the estimation procedure to the lack of uniqueness of the copula given that ∆ is a discrete (binary) variable. In this situation, to quote Embrechts (2009), "everything that can go wrong, will go wrong". Details about copulas for discrete variables may be found in Genest and Nešlehová (2007).
Instead, the idea is to work on the joint conditional expectation so as to bypass the issues related to the copula of (Y, ∆, X). In short, our intention is to discard the problem by obtaining the copula of (Y, X) conditionally on ∆ = 1, for which no specific technical difficulties are involved. To that end, using the notations H u (resp. h u ) as a shorthand for a given distribution (resp. density) conditionally on ∆ = 1, first note that Y X and f u denote the conditional densities of (Y, X) and X given ∆ = 1, respectively. Hence, using the definition of a copula function in a similar spirit as Noh et al., one obtains Applying this equality in the context of quantile regression, interestingly, one eventually retrieves an expression analogous to (2.2): ). Note that the copula in question in (2.8) is determined by strictly fully observed data. Hence, standard literature on copulas can be manipulated without any censoring related constraints. Given estimators p . . , d, satisfying certain high-level conditions which will be given in Section 3, this suggests to estimate the quantile regression in the presence of censoring by the empirical analogue of (2.8), that is , and where p c u Y X denotes an estimator of c u Y X based on the four-step procedure described in Section 2.1. Explicitly, where any two-dimensional kernel density estimator may be used for each bivariate copula density p c u Y X j , j = 1, . . . , d, such that condition (C7) of Section 3 holds, and where p c u X 1 ...X d |Y is estimated by standard parametric vine procedures.
The resulting quantile regression estimator in (2.9) may then be viewed as a simple weighted quantile of the observed response variable, and is therefore easy to implement in practice using the efficient quantile regression code developed by Portnoy and Koenker (1997) and Koenker (2005). Nonetheless, in the context of multivariate covariates, the estimation of G C (·|x) requires further assumptions to overcome dimension related issues. Popular choices in the literature include, among others, independence between C and X, the Cox model or the single-index model on C|X = x. Illustrations of such assumptions are treated in our simulation study.
As an interesting property, and similarly to the case without censoring, we note that the obtained regression function estimator is automatically monotonic accross quantile levels. Applying analogous arguments to the ones adopted in the proof of Theorem 2.5 of Koenker (2005), one can indeed determine that Conclusively, in parallel to what has been stated for the uncensored case, the resulting estimator p m τ (x) defines a rich class of estimators built on the many different existing methods available in the literature for estimating copula densities and marginal distributions of both complete and censored data.

Asymptotic Properties
We establish in this section the asymptotic normality of the proposed estimator p m τ (x). To that end, we first report the set of regularity conditions as well as the required high-level conditions on all estimators involved in the expression of p m τ (x). We then develop an asymptotic representation of our estimator for a general d-variate covariate. As the latter will result in a somewhat unpleasant expression for the asymptotic bias and variance for a general multivariate covariate, and given that the analytical reasoning is similar in spirit, we eventually restrict ourselves to the detailed asymptotic expression for the case d = 2.
For a fixed but arbitrary point of interest x in the support of X, denoted by supp(X), let us suppose that there exists a neighborhood V F u (x) of F u (x) such that the following regularity conditions hold: (C1) The conditional distribution F T |X of T given X admits a conditional density f T |X that is continuous, strictly positive and bounded uniformly on R × supp(X).
Concerning the high-level conditions, it is assumed that the multivariate copula density c u Y X is estimated using the proposed four-step strategy of Section 2.1, and that, for simplicity, the d bivariate kernel copula estimators of step (2) are based on the same bandwidth H = h 2 I for a certain h > 0. The following conditions are then assumed to hold: (C5) The marginal c.d.f. estimators are such that: and τ G C = sup{t : G C (t) < 1}.
(C7) The multivariate copula estimator is such that: where ∂ j denotes the partial derivative with respect to the j-th argument.
Assumption (C1) is standard in the context of quantile regression estimation. As for condition (C2), this is similar to assumption (C3)-(i) in Noh et al. (2015) for the simplified case with no censoring, with an additional requirement on the conditional censoring probability that is resulting from the initial transformation of synthetic observations. Assumption (C3) is likewise emanating from the handling of censoring through these observations, and is rather usual in survival analysis. Note that, in the quantile regression framework, the latter assumption amounts to defining a natural upper bound for the quantile of interest that can be studied. Assumption (C4) reports a set of technical conditions to be met.
As regards conditions (C5)-(C7), (C5) is routinely made in the copula framework. For instance, it is readily satisfied for the empirical distributions when only uncensored observations are taken into account, and their rescaled versions which are prominent in the copula literature. Assumption (C6) imposes restrictions on the estimator one may consider for the conditional distribution of the censoring variable and is, for instance, fulfilled for a simple Kaplan-Meier estimator for G C (see e.g. Theorem 2.1 in Chen and Lo (1997) for sufficient and necessary conditions for (C6)). Lastly, the uniform consistency of the kernel density estimator required by assumption (C7) is, for instance, alluded to in Geenens et al. (2014) for the probit-transformed copula estimator.
We now state the main result of this section that holds for a general d-dimensional covariate vector and for all bivariate kernel copula estimators based on the same bandwidth h. In practice, however, it may be recommended to adopt an unconstrained and non-diagonal bandwidth matrix, as is detailed in Section 4 of Geenens et al.. Nevertheless, when considering this general situation, the theoretical results become less tractable while equivalent in nature to the simplified situation considered here.
Theorem 3.1. Let h ≡ h n → 0 be the common bandwidth of the d bivariate kernel copula density estimators. For h satisfying nh 2 → ∞ as n → ∞, and under assumptions (C1)-(C7), we havè Theorem 3.1 implies, quite naturally, that the asymptotic behaviour ofm τ (x) will be characterized by the properties of the copula estimator, specifically through its nonparametric feature, provided that the estimation of p G C (·|x) is 'reasonable' when confronted to a multidimensional covariate vector (assumption (C6)). In particular, this suggests that the detailed discussion of Geenens et al. about the asymptotic bias and variance of their distinctive bivariate copula estimators may be transcribed in our context. Additionally, Theorem 3.1 also covers an asymptotic representation of the copula-based quantile regression estimator when all responses are fully observed. In this situation, one would indeed obtain a similar result for the proposed semiparametric procedure, with the removal of all censoring related terms, that is w(x), W i (x), i = 1, . . . , n, and the superfluous conditioning on ∆ = 1 for the copula densities and marginal distributions.
We now consider a detailed asymptotic representation of our estimator for the simplified case where d = 2, and for a general nonparametric estimator of the bivariate copula densities. For convenience, we use the notation c u k as a shorthand for c u Y X k , and similarly for other functions depending on (Y, ∆, X k ), k = 1, 2.
Corollary 3.2. Suppose that the assumptions of Theorem 3.1 hold for the case d = 2. Furthermore, suppose that the bivariate nonparametric copula estimators of c u 1 and c u 2 are such that Suppose furthermore that the following technical conditions hold:

Then, the copula-based quantile regression estimator at any point of interest
Corollary 3.2 reports the asymptotic normality of our estimator at the expected convergence rate, implied by the nonparametric estimation of the bivariate copula densities. Depending on the choice in step (2) of the kernel density estimator fulfilling the conditions of Corollary 3.2, simple plug-in of the expression of Z nj k , k = 1, 2, in all quantities built upon the latter may then lead to the detailed, although arduous, expressions of the asymptotic bias and variance of the proposed estimator.
Furthermore, as this had yet to be covered, it is worth stressing out that Corollary 3.2 also encompasses the asymptotic normality of the suggested estimator based on semiparametric vine copulas with strictly complete data. Similarly to what has been stated for Theorem 3.1, one is indeed only required to withdraw all censoring related elements from Corollary 3.2 to obtain the expressions of the asymptotic bias and variance of the proposed semiparametric quantile regression estimator for complete observations.

Simulation Study
In this section, we assess the practical finite-sample performance of the proposed methodology by means of Monte Carlo simulations. For this purpose, we first present a brief numerical study to further motivate the semiparametric copula strategy we intend to adopt for multivariate problems for both complete and censored observations. Secondly, focussing on survival data, we illustrate the flexibility of our estimator based on the proposed copula modelling by showing promising results with respect to competitors, including when the generated scenario is to the advantage of the latter. All the simulations are carried out using the statistical computing environment R (R Core Team (2014)) and its freely accessible packages.

Assessing the semiparametric copula estimation
This first section aims at numerically illustrating the choice of our semiparametric copula estimation strategy. For this purpose, we consider two distinctive data generating processes (DGP) and compare our methodology with fully parametric and nonparametric procedures one might consider for the estimation of a multivariate copula density. For the general simulation settings, we consider B = 500 repetitions of each DGP; three (average) levels of censoring (0%, 30% and 50%), three sample sizes (n ∈ {100, 200, 400}) and the quantile level of interest τ = 0.3. As the object of interest here is the copula modelling, when censoring is introduced, we only consider the simple case of independence between the censoring variable and the covariate vector in order to keep the estimation of p G C needed for (2.9) uncomplicated, that is, using the Kaplan-Meier estimator. The detailed DGPs are as follows: • DGP A: (F T (T ), F 1 (X 1 ), F 2 (X 2 )) ∼ Gaussian copula with parameters (ρ T 1 , ρ T 2 , ρ 12 ) = (0.3, 0.9, 0.5). Given standard uniform marginal distributions for all three variables, the resulting true quantile regression may be calculated as m τ (x) = Φp −0.2Φ −1 (x 1 ) + Φ −1 (x 2 ) + 0.4Φ −1 (τ )q (see Noh et al. (2015)). To include censoring, we introduce the variable C ∼ U [0, M ], where the parameter M is computed in order to obtain the desired average censoring proportion (M = 5/3 for 30% and M = 1 for 50%).
For any general copula-based regression estimator, the marginal distribution estimations are performed, as suggested in Section 2.1, using rescaled versions of the empirical distributions: where n u = n i=1 ∆ i is the number of uncensored observations. For the distinctive copula-based estimators, we consider the following procedures: p m SP cop, τ : semiparametric estimation strategy detailed in Section 2.1. That is, we first estimate the d bivariate copulas of interest employing the probit transformation technique of Geenens et al. (2014) coupled with local likelihood estimation based on quadratic polynomials. To that end, we follow their proposed nearest-neighbor bandwidth selection procedure. Concerning the estimation of the d-dimensional 'noisy' copula density (2.5), as mentioned above, we apply standard vine techniques using the R package VineCopula. Specifically, we adopt one automatically selected tree structure for the simplified decompositon of the copula density among many R-vine candidate structures (see Dißmann et al. (2013)), and subsequently determine the appropriate pair-copula family to be selected and parametrically estimated. The selection criterion for bivariate copulas is chosen to be the Akaike information criterion (AIC), which revealed to be adequate in the R-vine context (see Brechmann (2010), chap. 5), and ten potential family candidates, together with their rotations, are considered: eight of them are Archimedian (Clayton, Gumbel, Frank, Joe, Clayton-Gumbel, Joe-Gumbel, Joe-Clayton and Joe-Frank), and the last two are elliptical (Gaussian and Student t).
p m N P cop, τ : fully nonparametric estimation of the d-dimensional copula using vine techniques. Specifically, while the vine structure is kept identical, here all bivariate building blocks are estimated using the local likelihood technique based on probit-projected data with, here again, the bandwidth selection procedure of Geenens et al. (as is studied in Nagler and Czado (2015)). Given its fully nonparametric nature, it should be mentioned that this estimator is not covered by the theoretical results of Section 3.
p m P cop, τ : fully parametric estimation of the d-dimensional copula density, where all bivariate copulas are estimated using the previously mentioned candidate families and selection criteria. However, unlike the above-mentioned estimators, we do not force here any structure for the vine decomposition. As a consequence, no explicit distinction is imposed between dependence of interest and noisy dependence. Instead, one data-driven selected structure is adopted, regardless of the arguments of Section 2.1. This will allow us to analyse the impact of such dependence distinction in our regression context, as is discussed below. Finally, as it is the case for p m N P cop, τ , this estimator is not covered by the asymptotic theory of Section 3.
Both DGPs concentrate on the situation where the dependence structure between the response variable and the covariate vector is characterized by a parametric copula. In such circumstances, p m P cop, τ will have a critical advantage, and may serve in order to evaluate the impact of the nonparametric part of the estimation scheme, especially when the dimension of the covariate vector increases. As a performance criterion, we consider here the empirical integrated mean squared error (IMSE), defined as where {x i , i = 1, . . . , N } is a generated random sample of size N = 10 serving as an evaluation set spread on the domain of X, and p m (b) τ (·) denotes the regression estimation for the b-th simulated sample.
The results of our simulation study are summarized in Table 1. Based on these, we detail our analysis in two parts, as the outcomes of our study offer relevant information on both the copula decomposition choice and the type of bivariate estimators one may adopt in the multivariate setting. Note that, for both DGPs, as the dependence structure is specified by a Gaussian copula, the simplifying assumption intrinsic to the vine decomposition is here applicable (see Theorem 4 in Stöber et al. (2013)). In our context, this means that any observed difference between copula strategies is not to be attributed to a possible violation of the underlying simplifying assumption.
Focussing first on our decomposition strategy, we note that, as expected, p m P cop, τ globally outperforms p m SP cop, τ and p m N P cop, τ . However, strikingly enough, this is not observed for DGP A for different censoring proportions, where p m SP cop, τ details better results. This is interpreted here as evidence for the validity of our arguments regarding the decomposition choice: as the censoring proportion grows, the number of observations actually entering the copula estimation becomes more moderate, hereby implying two opposite effects in this context. First, the propagation of estimation approximations tends to be more important, signifying that the further we decompose, the more sensitive becomes the estimation of the involved bivariate copulas as these are tributary of the quality of previously estimated bivariate blocks. Using a purely data-driven decomposition may then result in a poor fit of the (conditional) copula of the response variable with one of the covariates, as it is not required that the latter would be primarily treated. This is interpreted as the reason why p m SP cop, τ is able to outperform p m P cop, τ , admittedly by a small amount, when censoring increases for a fixed sample size n ∈ {200, 400}, even though the simulated scenario is issued from a purely parametric copula. However, on the other hand, when observations are more scarse, it is well-known that nonparametric estimations become more sensitive than parametric counterparts. This explains why the estimation results for p m SP cop, τ are not superior to those of p m P cop, τ for n = 100 with 50% censoring, as the former requires the nonparametric estimation of two bivariate copulas, whose complexity compared to p m P cop, τ seems to override the positive effects of our decomposition choice. Overall, these noteworthy results for DGP A illustrate the effectiveness of our proposed copula decomposition in the regression context. When augmenting the covariate vector dimension, the price of estimating now three nonparametric bivariate copulas quite logically exceeds the potential gain of concentrating efforts on the dependence of interest. This is identified in DGP B.
Concentrating now on the modelling choice for the noisy dependence, the comparison between p m SP cop, τ and p m N P cop, τ offers valuable information, particularly in DGP A, as the only distinctness here is the estimation of a unique bivariate copula density c u X 1 X 2 |Y . Visibly, the implication of keeping a nonparametric approach for this part seems to be rather severe. This finding clearly also applies when the dimension of the covariate grows.
Conclusively, the recommended copula modelling strategy seems to propose an adequate tradeoff between preventing the serious effects of a purely nonparametric estimation and providing the flexibility needed to overcome possible shortcomings associated to a purely parametric alternative.

Comparison with other estimation methods
The objective of this second section is to provide a comparison study between the proposed copulabased methodology and existing competitors for survival data. On a general note, given that for multidimensional covariates an appropriate estimation of the conditional distribution G C (·|x) for the weights W (x) in (2.9) may be of crucial influence, we consider distinctive scenarios contrasting in the impact on p G C (·|x) in order to provide a sufficiently broad view on the performance of our methodology.
We examine 2 general simulation models, with B = 500 repetitions of each; two (average) levels of censoring (30% and 50%), two sample sizes (n = 200 and n = 400) and four values for the quantile level of interest (τ ∈ {0.1, 0.3, 0.5, 0.7}). Specifically, we consider general data arising from a Cox regression model. Based on the hazard function (defined as h(t) = f T (t)/(1 − F T (t))), this prominent model for analysing survival times specifies that h(t|X i ) = h 0 (t) exp(β T X i ), i = 1, . . . , n, where h 0 (t) is the so-called baseline hazard function. In this instance, given a nonnegative time-to-event variable, simple algebraic manipulations show that the general conditional quantile regression can be written as m τ ( For every proposed setting, we choose β = (1, −3/4, 1/2, 1/4, −3/5) T and d = 5. Furthermore, the baseline distribution is set to be standard exponential, and, consequently, for a given vector x the τ −th conditional quantile function is given by The distinction between our simulation scenarios is to be found in the covariate vectors that are taken into account. Let us consider 5 covariates X j , j = 1, . . . , 5, with standard uniform distributions. In order to allow for covariate dependence, we simulate (X 1 , . . . , X 5 ) form a 5-variate Gaussian copula, with parameters (ρ 12 , ρ 13 , . . . , ρ 45 ) = (0.3, 0.4, 0.5, 0.6, 0.7, 0.3, 0.4, 0.5, 0.6, 0.7). To distinguish our scenarios, we consider two covariate vectors given by X (1) = (X 1 , X 2 , . . . , X 5 ) and X (2) = (X 1 , exp(X 2 ), X 3 , X 4 , X 5 ), for which the resulting quantile regressions are both singleindex models in X (1) and X (2) , respectively. The motivation of our simulation study is to highlight the performance of our procedure when its competing estimators are both in an ideal situation and a slightly altered version of it.

Model 1: covariate vector X (1)
We first consider the simple case of data following a Cox regression model issued from covariate vector X (1) . The distributions and parameter values of the censoring variables define the following scenarios: • DGP C: C ∼ Exp(λ C ), with λ C independent of X (1) . To attain the desired average censoring proportions, we fix λ C = 0.464 and λ C = 1.083 (corresponding to approximately 30% and 50% censoring in average, respectively). The exact conditional censoring probability given X (1) = x is calculated as λ C /(λ C + exp(β T x)) and is hence a decreasing function of β T x, making us conjecture better results for higher values of β T x. Note that in this scenario, G C (·|x) boils down to G C (·), for which the Kaplan-Meier estimator is suited.
• DGP D: C ∼ Exp(λ C (x)), with λ C (x) = 3/7 × exp(β T x) and λ C (x) = exp(β T x) corresponding to 30% and 50% censoring, respectively. In this scenario, the conditional censoring probability is independent of x, and an adequate estimation of G C (·|x) is fulfilled with a Cox model hypothesis for the relationship between the covariates and the censoring variable itself.
For comparison purposes, we consider the four following estimators: p m cox, τ : parametric estimator exploiting the information related to the parametric Cox model setting. This estimator will serve as a reference for the ideal, yet unknown in practice, situation.
where p β is estimated by maximum partial likelihood, and p λ T 0 is the maximum likelihood estimator of the exponential baseline distribution.
p m si, τ : Single-index regression estimator studied in Bücher et al. (2014), where the censoring distribution is supposed to be independent of x. For the univariate nonparametric part of the estimation process, 10 different bandwidths are selected (h ∈ {0.05, 0.1, . . . , 0.5}), and the optimal choice is performed using the described leave-one-out cross-validation procedure.
Note that the quantile regression of interest m τ (x) here is indeed a single-index model in (X 1 , X 2 , . . . , X 5 ). This should provide a critical advantage to the performance of p m si, τ .
p m cop, τ : our copula-based estimator from estimation equation (2.9), where the conditional distribution G C (·|x) is supposed to be independent of x and thereby estimated by the classical (unconditional) Kaplan-Meier estimator.
p m (cox) cop, τ : our copula-based estimator, where the relationship between the covariates and the censoring variable is supposed to follow a Cox regression model. Namely, p G C (c|x) is estimated by Following the arguments of Section 2.1 and the results of Section 4.1, both copula-based procedures are implemented employing a semiparametric estimation for the d-variate copula built on the aforementioned candidate families and selection criterions.
In order to compare the studied estimators' performance, we consider an integrated version of the median absolute estimation error, that is where p m τ is a generic estimator of m τ , {x i , i = 1, . . . , N } is an evaluation set corresponding to a generated random sample of size N = 10, spread on the domain of X (1) , and med (B) denotes the median taken over all B = 500 simulations. The choice for this robust L 1 -type of measure is motivated by the fact that the optimization routines involved in the single-index procedure may yield very unlikely results with a small probability, hereby strongly disadvantaging the estimator when considering a L 2 -type of error measure. The same reasoning is underlying the determination of a robust dispersion measure on estimation errors, taken as the averaged interquartile range. More precisely, we choose and Q (B) 1 stand, respectively, for the third and first quartiles taken over all simulations. The results of our simulation study for this model are reported in Tables 2 and 3 for 30% and 50% of censoring, respectively. As expected, the Cox regression estimator p m cox, τ , serving as a reference case here, outperforms the other estimators for both scenarios and every sample size, quantile level and censoring percentage. More interestingly, while the single-index estimator p m si, τ quite logically displays better overall results than our copula-based estimators, it is worth noticing that the difference seems to fade away when moving to higher quantile levels, especially for higher levels of censoring percentage. This is considered as a first encouraging result for the copula-based estimators as the simulated model is very strongly to the advantage of p m si, τ . This indicates that our proposed procedure is flexible enough to compete with p m si, τ in its ideal setup when the number of observations actually entering the estimation scheme becomes moderate, that is when censoring percentage is important and the quantile level of interest is high. Of course, if this is not the case, there is no a priori reason to believe that the copula-based estimators could outperform p m si, τ when the latter is considered in its optimal setting. Figure 2 serves to illustrate these results for the particular case of DGP C with n = 400 and for high quantile levels (τ ∈ {0.5, 0.7}). For the sake of brevity, we only report here a graphical comparison between the copula-based estimators and the single-index estimator, as p m cox, τ clearly outperforms the latter. As can be observed from Figure 2, the difference between p m si, τ and both copula-based estimators is graphically modest, especially for the highest quantile of interest τ = 0.7, hereby reinforcing the previously discussed results of tables 2 and 3. Furthermore, note that all estimators tend to present worse performances for low levels of β T x. This is expected in this simulation setting, as these levels correspond to the region of β T x where the exact conditional censoring probability is the highest.
Focusing again on tables 2 and 3, we further note that, while there is logically a difference between the performances of p m cop, τ and p m (cox) cop, τ depending on the simulated scenario, the effect of an appropriate modelling of p G C (·|x) seems to be rather limited in this simulation setup for our estimators. Additionally, while both p m si, τ and p m cop, τ rely here on the assumption of independence between the censoring variable and the covariate vector, the latter seems to numerically behave better when confronted to the violation of this hypothesis. This can be observed from the comparison between DGP C and D of the dispersion measures of both estimators, once again especially for high levels of quantile values and censoring percentages. Of course, one could argue that, given the results for the dispersion measure of p m si, τ for high quantiles, this can partly be due to a poor smoothing parameter choice. However, the latter was implemented using the proposed methodology of Bücher et al. (2014), just as a practitioner would have resolved to act.

Model 2: covariate vector X (2)
To complete our simulation study, we now consider in this section the second model where the time-to-event data is simulated using covariate vector X (2) . Consequently, the resulting quantile regression is no longer a single-index nor a Cox regression model in (X 1 , X 2 , . . . , X 5 ). For the sake of brevity, we only consider here the situation where the censoring variable is independent from the covariate vector, as the impact of a dependent scheme has been treated in model 1. Specifically, we simulate the censoring variable from an exponential distribution with parameter values 0.208 and 0.486 for approximately 30% and 50% censoring. For comparison purposes, we consider the four estimation procedures described in model 1, given covariate vector (X 1 , X 2 , . . . , X 5 ). The copulabased estimation procedures are constructed using the same semiparametric modelling strategy as for model 1.
The resulting performance of the considered estimators are depicted in Table 4, once again in terms of IM AE. In this simulation scheme, we observe that both copula-based estimators tend to outperform their competitors, with the exception of estimations for very low quantiles of interest, where the effect of a ('small') misspecification of the underlying model seems to be moderate for both p m cox, τ and p m si, τ . As we move away from these low quantile levels, the consequences of misspecifying the model become more severe for the competing estimators, notably for a purely parametric approach ( p m cox, τ ). In contrast, the copula-based approach presents satisfactory results for varying censoring proportions and sample sizes, especially when keeping in mind that the simulated scenario is 'only' a slightly altered version of the ideal scenario for its competitors. Conclusively, this advocates, here again, for the flexibility of our procedure and its withstanding to misspecification of the underlying model for multidimensional problems when comparing with other semiparametric or fully parametric modelling techniques.

Real Data Application
We present in this section a brief application of our procedure by analysing the Colorado Plateau uranium miners cohort data (see e.g. Lubin et al. (1995), Langholz and Goldstein (1996)). The object of the study, for which 3347 Caucasian male miners having worked at least a month in the uranium mines of the Colorado Plateau were followed, is to investigate the risk of lung cancer related to smoking and radon exposure. Hence, the event of interest is defined as the time till lung cancer death (expressed as the logarithm of number of years), which affected a total of 258 miners. Besides failure time, the study also includes information about age at entry to the study, cumulative smoking (in number of packs) and radon exposure (in working level month (WLM)). As the original data set is prone to heavy censoring (92.3%), and given the illustrative nature of this section, we first define a subsample on which the analysis will be performed. To that end, in order to preserve as best as possible the nature of the population at risk, we decide to define a threshold on the radon exposure above which observations will enter the subsample. The value of the threshold is practically chosen as a trade-off between censoring proportion and actual number of observations that are to be selected. Specifically, we find that by defining a threshold of 2831 WLM on radon exposure, a subsample of 176 observations is constituted, 55 of which were subject to the event of interest. Scatterplots of the selected data are represented in Figure 3, where X 1 is the age at entry into the study, X 2 is the cumulative radon exposure and X 3 is the cumulative smoking.
Regarding the data analysis, endorsing the role of a practitioner, we are faced with the choice of an appropriate estimator for the application of quantile regression in this context. We therefore consider the following distinctive candidates: where p β is estimated by maximum partial likelihood, and p H 0 is the Nelson-Aalen-type estimator of the cumulative baseline hazard.
As a general evaluation measure to compare models for the present data set, we consider the median quantile loss from predicting the τ -th conditional quantile of T for the uncensored observations. In other words, we use the following cross-validated prediction error criterion: where p m −i τ denotes any estimator of m τ based on all observations except the i-th one. For the implementation of the copula-based regression estimator, we adopt the methodology of Section 2.1 and our simulation study by opting for a semiparametric modelling of the fourvariate copula density, where the bivariate copulas of the response variable with each covariate are estimated using the procedure of Geenens et al. (2014) with quadratic polynomials along with the proposed data-driven bandwidth selection scheme. The remaining trivariate noisy copula is, afterwards, estimated using vine techniques with the same candidate families and selection criterion as in Section 4. Additionally, a general appropriate hypothesis has to be made for the modelling of the conditional distribution G C (·|·). As it is shown in the literature that independence is unsuitable for the data of the study (see e.g. Leng and Tong (2013)), we propose to model the conditional distribution through a semiparametric Cox regression for the censoring time with respect to the covariates. Therefore, as is common in practice, the parameter of the regression is estimated by maximum partial likelihood while the Breslow estimator is used for computing the baseline survivor function.
Concerning p m si,τ , the bandwidth choice is performed using the proposed leave-one-out crossvalidation procedure on normalized covariates with 15 candidates (h ∈ {0.1, 0.2, . . . , 1.5}). Note that, as the latter procedure is computationaly costly, it is here assumed that the selected bandwidth is constant for all datasets entering the estimation of PE( p m si,τ ), which is of little impact in this context given that only one observation at the time is to be removed from each dataset.
The results of the evaluation measure are reported in Table 5 for four quantile levels of interest τ ∈ {0.1, 0.3, 0.5, 0.7}. It is observed that, in terms of cross-validated prediction error, our copulabased approach depicts quite confidently the best performance for every considered quantile level of interest. By comparison, the single-index structure seems to be inappropriate as such for the studied data set, and should possibly require further work and attention on, for instance, transformations of covariates. Furthermore, its sensitivity to higher quantile levels is, here again, highlighted. Lastly, despite being part of every practitioner's toolbox for survival analysis, the proportional hazards regression estimator is also relatively largely outperformed by the increased flexibility of the copula-based approach. Hence, this example clearly illustrates the ability of our estimator to adapt to the underlying regression structure of the data.  Conclusively, recalling that, from equation (2.10), our procedure enjoys the additional valuable property of being automatically monotonic across quantile levels, this real data application plainly highlights the relevance of our estimator for flexible analyses of multivariate censored data.

Conclusion
In this work we have proposed a semiparametric copula-based quantile regression estimator in the context of potentially right-censored responses. On a general note, for data with or without censoring, and motivated by the regression context, a novel semiparametric estimation approach for the implied copula density was studied. Furthermore, in parallel to the procedure of Noh et al. (2015), the proposed regression estimator in this work is obtained as a weighted quantile of the observed response variable, hereby opening the door to the practical use of the quantile regression code developed by Portnoy and Koenker (1997) and Koenker (2005). Asymptotic normality of the resulting estimator for both complete and censored data was obtained with convergence rate determined by the nonparametric estimation of bivariate copula densities. Finally, supporting the objective of this work, an extensive simulation study and a real data application have been carried out to illustrate both the validity of the semiparametric copula estimation and the increased flexibilty of our procedure in comparison with existing alternatives, especially when no a priori knowledge about the functional form of the quantile regression is available in practice.

Appendix
We develop in this appendix the proofs of Theorem 3.1 and Corollary 3.2.

Proof of Theorem 3.1
Define , and a n = ? nh 2 . Observe first that, by definition of p m τ (x), a n (m τ (x) − m τ (x)) = arg min s p A n (s).
Furthermore, given that ρ τ is a convex function and that x W i (x) p c u Y X p p F u Y (Y i ), p F u (x)q ≥ 0 for all i = 1, . . . , n, we have that s → p A n (s) is convex. The idea is then to develop an expression of p A n (s) leading to the application of the quadratic approximation Lemma of convex functions (Basic Corollary in Hjort and Pollard (1993)). To that end, we have that (Knight's (1998) and 0 ≤ R(u, v) ≤ |v|. Hence, we may write p A n (s) = −s p A 1n +Â 2n (s), Focusing first on p A 2n (s), we will show that where f T |X is the conditional density of T given X and w(x) = p(x)/rP(∆ = 1)c u X (F u (x))s. To that end, we write p A 2n (s) = A 2n (s) + ∆ 1n (s) + ∆ 2n (s), Concentrating on each term, we will show that For the proof of (A.2), we first establish that E(A 2n (s)) = n E pR( 1 , s/a n ) where F T |X denotes, as in Section 2.1, the conditional c.d.f. of T given X, and where, for the second equality, we used (see Section 2.2) the fact that, for any measurable function ϕ : R → R, for some θ ∈ (0, 1), with R(t) = f T |X (m τ (x) + θt|x) − f T |X (m τ (x)|x). This yields To that end, observe that, for n sufficiently large, and under assumptions (C1),(C2) and (C3), , as a n /n converges to 0.
We will show that all these quantities converge to 0 in probability. Starting with ∆ 3n , we have, for a large n, under the condition that assumptions (C2), (C4)-(i), (C6) and (C7)-(i) are met.

Proof of Corollary 3.2
Given that, by the result of Theorem 3.1, the asymptotic behavior of p m τ (x) will be dictated by the expression of the multivariate copula estimator, we will concentrate on the latter. Using the asymptotic expressions of both bivariate copulas along with the multivariate copula decomposition