Traditional regression models focus on estimating the conditional mean. However, relying solely on the mean is inadequate since real data often exhibits heteroscedasticity. Instead, there is a growing recognition of the need to comprehend the underlying distribution. Ideally, estimating the entire conditional distribution would be beneficial. Foresi and Peracchi (1995) and Klein et al. (2015) suggested distributional regression methods to capture the full conditional distribution. However, they can be computationally demanding and pose challenges in terms of interpretability.

Besides, quantile regression (QR) (Koenker and Bassett, 1978) provides an efficient method to achieve a comprehensive understanding of the data. It is particularly robust in handling data with high skewness, non-Gaussian errors, and outliers. We can explore the conditional distribution through QR without making assumptions regarding the global distribution (Chapter 2.5 in Koenker, 2005). Many fields that require knowledge beyond the conditional mean utilize QR these days. For example, meteorology adopts QR to forecast heavy rainfall, healthcare examines treatment side effects (Powell, 2020), and finance calculates the value-at-risk of a risky investment (Keilbar and Wang, 2022) using QR.

Estimating quantiles for complex data requires nonlinear models. Tree-based models targeting quantiles include quantile regression forest (QRF) (Meinshausen and Ridgeway, 2006) and quantile extremely randomized trees (QERT) (Geurts et al., 2006). Boosting models such as gradient boosting machines (GBM) (Friedman, 2001), LightGBM (Ke et al., 2017), and extreme gradient boosting (XGBoost) allow for quantile estimation by substituting their original loss functions with a check loss. Quantile smoothing splines (QSS) (Koenker and Portnoy, 1994) and support vector quantile regression (SVQR) (Takeuchi et al., 2006; Li et al., 2007) are another variants of nonlinear quantile regression models.

The crossing of the fitted quantiles is often an issue in the quantile regression models. To prevent the fitted quantiles from crossings, researchers often impose constraints (He, 1997; Liu and Wu, 2011; Shin and Jung, 2023). Although the proposed model does not employ certain constraint, we observe that crossings can be much relieved as a byproduct.

In this paper, we aim to apply the concept of composite quantile regression (CQR) (Zou and Yuan, 2008) to various ML models. CQR combines objective functions for multiple quantiles and minimizes them simultaneously. Bloznelis et al. (2019) and Zhao and Xiao (2014) showed that the combined estimator has a higher efficiency than the original quantile regression. However, CQR primarily treat the central tendency of linear models, which does not well matched for the terminology ‘quantile’. Xie and Wen (2019) and Hatalis et al. (2019) integrate CQR concepts with deep learning models, but those algorithms still primarily target the mean or median. One of the assumptions of CQR is that all quantiles share the same regression parameter, and only varying the intercept term. This assumption is restricted when analyzing heterogeneous data. Therefore, our goal is to relieve this assumption and to develop methods that can estimate all range of quantiles for nonlinear and heterogeneous data.

We propose a new methodology called local quantile ensemble (LQE). LQE aims to ensemble models for several quantiles in a straightforward way, averaging multiple quantile fits near a target. This method is not restricted to specific algorithms, making it applicable to any quantile-targeted models. While traditional quantile estimation focus on single quantile, LQE averages multiple quantiles near the target. As a result, its performance may improve through the bias-variance trade-off. Additionally, LQE can reduce quantile crossing, which is one of the main concerns in QR.

Our paper will proceed as follows: First, we review some quantile machine learning algorithms in Section 2.1. In Section 2.2, we introduce the details of LQE, and illustrate the bias-variance trade-off with a pilot example. Sections 3 and 4 show extensive numerical and visual results of the simulated and real data sets, respectively. Monte Carlo simulations and real data analysis reveal the usefulness of the proposed methods.

2.1. Review of quantile machine learning methods

Consider a univariate response variable y ∈ ℝ, and let x represent covariates with p dimensions. Suppose we have independent samples (x_i, y_i), i = 1, . . . , n. The goal is to estimate the conditional quantile of y given x. Denote the true quantile as Q_τ(y|x), where τ ∈ (0, 1) indicates the quantile of interest. Then, true conditional quantile is defined as

Here, F denotes a conditional distribution of y given x. To obtain the estimates of (2.1), Q̂_τ(y|x), there are several classes of machine learning models.

Boosting models (e.g., GBM, LightGBM, and XGBoost) have been used mainly for mean regression, but they can be easily implemented for quantile estimation. Conventional boosting models employ mean squared error (MSE) or mean absolute error as a loss function. We substitute these with the mean check loss function to estimate conditional quantile. The check loss is defined as L_τ(y, ŷ) = |y – ŷ|{τ – I(y < ŷ)}, where I is an indicator function that returns 1 if the condition inside is true or 0 otherwise.

With the mentioned loss function, boosting models learn by gradually integrating weak learners. The weak learner is a model that performs relatively poor. Its performance is just above the random guess. Denote the number of iterations as M. Then, we can formulate boosting models for quantiles as follows:

The λ represents a learning rate, which is typically a very small constant. ĝ^m(x) and y_i,(m) are the decision tree and the i^th residual in the m^th iteration, respectively, for m = 1, . . . , M. Each tree is fitted to correct the previous errors. The algorithm stops either when the maximum number of iterations is reached or when a new decision tree no longer improves the loss function.

Meanwhile, QRF and QERT estimate quantiles by obtaining an empirical conditional distribution. These models calculate the weighted average of I(y_i ≤ y). Suppose we construct a QRF with H trees, where H is a positive integer. Following the notation of Meinshausen and Ridgeway (2006), let θ_h be a random parameter vector of the h^th tree for h = 1, . . . , H. Let R_l(x,θh) be the rectangular subspace generated by a single leaf l(x, θ_h). Then QRF combines all trees by averaging the weight vectors from h = 1, . . . , H trees as follows:

Note that the weight vector w_i(x, θ_h) indicates a proportion of observations that are contained in the leaf for x. We can subsequently obtain the estimates Q̂_τ(y|x) by substituting the estimates F̂(y|x) in place of F(y|x) in (2.1) as follows:

Once the empirical distribution of F̂ (y|x) is obtained, one can vary τ in (2.3) to target any quantiles of interest. Thus, we can obtain fits for multiple quantiles from a single fitted model. Furthermore, quantile crossing does not occur in QRF since the obtained empirical distribution is a monotonely increasing function of y.

QERT uses a similar fitting algorithm as QRF but amplifies the randomness. The primary differences are that each tree in QERT completely randomizes the cutpoints of the nodes and uses the entire sample to grow the tree. As a result, we can obtain the estimates using the same method as in (2.3), ensuring they do not exceed the above quantile fit. Further details of the algorithm can be found in Geurts et al. (2006).

Many studies compare various ML models to analyze real-world data with nonlinear patterns. To improve the performance ML quantile models, we suggest to combine multiple models with the ideas of CQR. The details of ensemble comes in the next section.

2.2. Local quantile ensemble

The methodology we propose here is local quantile ensemble (LQE). The core idea of LQE is to combine K estimates of quantiles 0 < τ₁< · · · < τ_K < 1 near the target quantile τ to yield a more accurate estimate. For simplicity, we refer to the K quantiles as local quantiles where the target quantile is located at the center. LQE combines estimates of the local quantiles around the target quantile as

We can easily implement LQE for models that return conditional quantile estimates, such as those in (2.2) and (2.3). While (2.4) can be extended to a weighted average, this paper only addresses the equally weighted version.

LQE has two advantages over previous studies. First, it is broadly applicable beyond linear models. Existing methods are largely based on the linear model, where the objective is to estimate the regression parameter. For example, Blozenlis et al. (2019) and Zhao and Xiao (2014) compute the average of the coefficient estimators for multiple quantiles. To extend this to ML models, we averages the fitted values instead of averaging the parameters. Second, LQE is capable of targeting any quantiles between 0 and 1. Most studies of CQR focus on estimating the central trend of data using full quantile information, typically limiting the target to the mean or median. In contrast, LQE can aims to estimate all conditional quantiles. As a result, LQE provides a more comprehensive approach to quantile prediction.

The ensemble effect can improve the performance of LQE. Ensemble approaches address problems where models may overfit or become trapped at a local optimum (Dietterich, 2000). Traditional quantile regression utilizes restricted information near the target quantile, making it sensitive to changes in observations above or below the fit (Koenker, 2005, Chapter 2.5). Consequently, estimates tends to be less accurate in regions with few samples or when τ is close to an extreme quantile (Li and Wang, 2019). By combining adjacent quantiles, LQE utilizes a broader range of samples, effectively addressing these issues. The bias-variance trade-off offered by LQE is crucial here. By averaging multiple quantile fits, LQE reduces the variance associated with individual model predictions, leading to more reliable and robust estimates. Despite an increase in bias, the overall performance of LQE often improves due to the substantial reduction in variance.

Next, we detail how to decide the local quantiles {τ₁, . . . , τ_K}. Two factors need to be specified: The range and the number of local quantiles. Combining quantiles that are far from the target quantile might not be helpful. Therefore, it may be necessary to restrict the range of the local quantiles. Denote the range as band(τ), which stands for bandwidth. The optimal length of band(τ) is influenced by the model and algorithm chosen, but establishing a theoretical optimum for each model is challenging. Therefore, we propose a simple heuristic: Use a wider bandwidth for quantiles with high densities. With the sparse density, each quantile fit is unstable. Even a small change in the target quantile can have a significant impact on the fit. Therefore, it is recommended to set band(τ) narrowly for those spaces. Because typical data sets have more samples near the center, we suggest the following rule for the bandwidth.

(Rule 1) The bandwidth for the target quantile τ, band(τ), is determined as

Note that the maximum length of the band(τ) in Rule 1 is 0.5, which is achieved when τ = 0.5.

Once we fix the bandwidth according to the target quantile τ, the next decision is the number of quantiles to combine, denoted by K. This is simpler than determining the range, as Koenker (2005, Chapter 5.5) notes that selecting a large K value does not significantly affect performance. Previous studies have often set K to 9 or 10 (Zhao and Xiao, 2014; Cannon, 2018). Furthermore, Zou and Yuan (2008) found no significant difference in the performance of the estimates beyond K = 9. When too many number of quantiles are averaged, the correlation between local quantile estimates increases. This hinders the variance reduction. Therefore, we fix K at 9 for LQE.

Pilot example: To demonstrate the bias-variance trade-off of the proposed method, we design a pilot experiment. Following Scenario 1 in Section 3, we generate the training data sets for R = 100 times. We fit a LightGBM model for target quantile τ = 0.1, 0.3, 0.5, 0.6, 0.8, and 0.9. We decompose MSE into the squared bias and variance with a fixed test data set of 1000 samples. Let Q^τ(r)(y∣xi) be the τ^th conditional quantile estimate given x_i in the r^th training data set. We calculate the squared bias and variance for the i^th test sample as follows:

Here, Q_τ(y|x_i) represents a true quantile, and Q̄_τ(y|x_i) is the average of the estimated quantile function over the R replication.

Table 1 summarizes the average of the estimated squared bias and variance across 1000 test samples. SQ refers to a single quantile from a model that fits the target quantile only. We observe that LQE reduces the variance while only slightly increasing the bias in general. The larger bias of LQE compared to SQ is a result of employing multiple quanitle fits. Nevertheless, our empirical findings indicate that the reduction in variance outweighs the increase in bias. This reduction in variance is particularly noticeable in the quantiles τ = 0.8 and 0.9, where the corresponding density of the error distribution is low. These results show the favorable effect of ensemble. In the next section, we will see that the variance reduction can lead to better performance in various ML models. Although we no more presents the decomposition of MSE, we believe the reduced MSE by the proposed methods is due to the bias-variance trade-off.

We compare the performance of our method to those of various ML quanitle regression methods. We consider two scenarios:

• (Scenario 1) We generate data following the Example 2 in Cannon (2018).

  y=(1-2x+4x2) exp (-2x2)+(0.2+0.08x) ɛ,

  where x is randomly drawn from U(−2, 2). We consider χ²(3) for the distribution of ε.

• (Scenario 2) We modify the Scenario 4 in Das and Ghosal (2018).

  y=sin(πx1)+cos(πx2)+(x12+1.5x22-0.5x2) ɛ,

  where the bivariate covariates x = [x₁, x₂] are randomly drawn from U(0, 1) × U(0, 1). The error term ε follows a t-distribution with 4 degrees of freedom.

Note that the generated data have a nonlinear trend and heteroscedastic errors in both scenarios. For the first scenario, we mainly focus on the error with the chi-square distribution, but we additionally consider the standard normal distribution for the graphical demonstration. For the second scenario, the original paper considered Gaussian errors, but we generate errors with a heavy-tailed distribution. The target quantiles are τ = 0.1, 0.3, 0.5, 0.6, 0.8, and 0.9.

Among ML algorithms capable of estimating conditional quantiles, we select GBM, XGBoost, LightGBM, QRF, and QERT. In addition, we fit QSS and SVQR for the visualization of Scenario 1 in Figure 1. We use R 3.5.1 program for QSS and SVQR, and Python 3.11 for the others. QSS is implementable with the rqss function in the package quantreg (Koenker et al., 2018). The kqr function from the package kernlab is applicable to SVQR. We can fit QRF and QERT using the functions RandomForestQuantileRegressor and ExtraTreesQuantileRegressor in the python package sklearn-quantile (Roebroek, 2022). We use the GradientBoostingRegressor function in scikit-learn (Pedregosa et al., 2011) for GBM, the lightgbm package (Ke et al., 2017) for LightGBM, and the xgboost package (Chen and Guestrin, 2016) for the XGBoost model. The LQE method can be easily implemented using existing packages.

We compare three models: Single quantile (SQ) and two versions of local quantile ensemble (LQE). The SQ is a genuine model with a single target quantile. LQE combines K = 9 local quantiles. To determine the magnitude of the local, we use Rule 1 in Section 2.2. Note that LQE requires to tune each of the K model. This increases the computing cost by K times. To reduce the tuning process, we consider a simpler version. First, we find the optimal parameter for the target quantile only. Then, we apply it to all the other local quantile models. The simpler one, LQE-S, needs a similar tuning time to SQ. Although the tuning time for LQE-S is much shorter than that of LQE, its results are quite comparable.

We split the generated data into training and test samples. The training data is used for model fitting and validation for parameter tuning, and the test data is used for measuring the prediction accuracy. In Scenario 1, the training and test sample sizes are 500 and 1000, respectively. In Scenario 2, we utilize a data set of 2000 samples for training and 5000 samples for the test. The simulations are repeated 100 times.

We use two criteria for evaluation. First, we calculate a mean square error (MSE) between the true quantile function and the predicted quantile fit over the test data, as follows:

This is a well-known indicator of the estimation accuracy. Next, we employ the crossing loss (Sangnier et al., 2016) to measure the undesirable quantile crossings. It is defined as

The crossing loss is nonzero when the quantile estimate for τ′ exceeds that for τ″ quantile. We record the crossing loss for only boosting models because QRF and QERT innately prevent quantile crossing. The crossing losses for these models are all zero.

We select the parameters through a grid search. The chosen parameters are the empirical minimizers of the mean check loss from 10-fold cross-validation. The range of parameters and their description are given in the appendix. We search for the same parameter space for SQ, LQE, and LQE-S within each ML model for a fair comparison.

Tables 2 and 3 show the results for Scenarios 1 and 2. Regardless of the ML model employed, LQE generally exhibits a lower MSE compared to SQ. It suggests that the ensemble effect of LQE results in improved performance. LQE-S, while slightly behind LQE in some cases, still often outperforms SQ, making it a viable option for scenarios where computational resources and time are constrained. The benefit of the proposed method is reassured in crossing losses. In most cases, LQE and/or LQE-S show much smaller crossing losses than SQ. In general, LQE provide a stabler fit which tends to keep the fits from crossing.

For graphical illustration, we apply quantile smoothing splines (QSS) and support vector quantile regression (SVQR) to Scenario 1. The algorithms for them are in Koenker and Portnoy (1994), Takeuchi et al. (2006), and Li et al. (2007). Figure 1 displays the fitted curves of SQ and LQE. Quantile crossings sometimes occur in SQ near the boundary, whereas they are rare in the proposed method. The numerical summary of Figure 1 is detailed in the appendix. Figure 2 illustrates the fitted curves of ML models in Scenario 2. LQE displays a slightly stabler surface than SQ.

In this section, we analyze three real data sets. We select data from three different fields. All data sets used here are available in the UCI machine learning repository. We briefly introduce each data set.

Bike: The bike-sharing data contains information on the number of bike rentals and the weather by day or hour gathered from the bike-sharing system. Among these, we use the data collected per day. Fanaee-T and Gama (2014) first analyzed this data, and Schallhorn et al. (2017) used D-vine quantile regression with it. The objective is to solve a regression problem predicting the total number of bikes rented per day. After excluding variables linearly dependent on the total number of rentals and the date variable, there are p = 11 variables and n = 731 observations. We randomly select 500 samples as training data and the others as test data.

Concrete: We analyze the concrete data to predict the compressive strength of the concrete based on its age and ingredients. The data was first proposed by Yeh (1998). Romano et al. (2019) conducted conformal inference using quantile regression. It contains n = 1030 samples, with p = 7 covariates. All covariates are continuous variables, such as the amount of fly ash and water in the concrete. We randomly divide the data into 700 training samples and 330 test samples.

Blog Feedback: The blog-feedback data was gathered by crawling and analyzing the raw HTML documents of blog posts. Buza (2013) initially collected the posts from 2010 to March 2012. Feldman et al. (2021) examined it to create prediction intervals using orthogonal quantile regression, and Alaa et al. (2023) analyzed it with conformalized unconditional quantile regression. We aim to predict the number of comments each post will receive during the next day. The data was already split into training and test sets, but we combine them. Final data has a dimension of n = 60, 021 and p = 280. Examples of covariates include the number of characters in a blog post and the date and time of publication. We randomly select 10,000 observations as training data and used the remaining observations as test data.

We fit the quantile ML models discussed in the previous section. The number of local quantiles, K, is 9 for LQE and LQE-S. All model fitting procedures are similar to those in Section 3. We search for the parameters in each model using a grid search. The details of the parameter spaces we search is in the appendix.

We measure the mean check loss, defined as

Tables 4, 5 and 6 summarize the check loss and crossing loss of the compared models from 100 random splits. The performance of LQE and LQE-S are frequently superior to that of SQ. These results demonstrate the proposed method can enhance the accuracy of predictions. LQE yields much lower crossing losses in most cases. LQE seems to suppress quantile crossing without additional constraints.

To demonstrate the computation efficiency of LQE-S, we also measure the computation time required for the parameter tuning process using the GBM model and the bike dataset. Because both SQ and LQE-S tune the target quantile only, their tuning times are roughly equivalent; thus, we only compare LQE and LQE-S. As shown in Table 7, LQE-S reduces the tuning time by approximately K = 9 times compared to LQE. This result highlights the substantial reduction in computing time achieved by LQE-S.

We introduce an ensemble learning method to improve quantile estimation and prediction of the existing machine learning (ML) models. LQE fits multiple local quantiles around a target and averages them. One can adapt LQE to a broad class of quantile regression models. It enhances performance by mitigating the large variance inherent in QR through the ensemble effect. Furthermore, LQE prevents overfitting and quantile crossing, as well as providing smoother fits.

However, LQE can be time-consuming due to tuning multiple parameters. To alleviate computing cost, we propose LQE-S, which tunes the parameter for the target quantile only and applies them to all other local quantiles. Our numerical analysis demonstrates that the practical performance of LQE and LQE-S are similar. In many cases, they achieve lower losses than single quantile (SQ) models. Moreover, the reduction in quantile crossing is evident. The wiggly fits and quantile crossings from the existing methods often significantly reduced.

This paper highlights the effect of ensembles near the target quantiles. We anticipate that LQE will perform effectively for another machine learning and deep learning models as well, which is an interesting area of our future research. In fact, the idea of LQE can be applied to the linear quantile regression models, and we developed early stage of theories.

Table A.1: The mean of MSEs at τ = 0.1, 0.3, 0.5, 0.6, 0.8, and 0.9 and crossing losses (with their standard error in the parenthesis) for 100 iterations under Scenario 1