Isotonic Distributional Regression

Abstract

1 Introduction

There is an emerging consensus in the transdisciplinary literature that regression analysis should be distributional, with Hothorn et al. (2014) arguing forcefully that

[t]he ultimate goal of regression analysis is to obtain information about the conditional distribution of a response given a set of explanatory variables.

Distributional regression marks a clear break from the classical view of regression, which has focused on estimating the conditional mean of the response variable in terms of one or more explanatory variable(s) or covariate(s). Later extensions have considered other functionals of the conditional distributions, such as quantiles or expectiles (Koenker, 2005; Newey & Powell, 1987; Schulze Waltrup et al., 2015). However, the reduction of a conditional distribution to a single-valued functional results in tremendous loss of information. Therefore, from the perspectives of both estimation and prediction, regression analysis ought to be distributional.

In the extant literature, both parametric and non-parametric approaches to distributional regression are available. Parametric approaches assume that the conditional distribution of the response is of a specific type (e.g. Gaussian) with an analytic relationship between the covariates and the distributional parameters. Key examples include statistically post-processed meteorological and hydrologic forecasts, as exemplified by Gneiting et al. (2005), Schefzik et al. (2013) and Vannitsem et al. (2018). In powerful semi-parametric variants, the conditional distributions remain parametric, but the influence of the covariates on the parameter values is modelled non-parametrically, for example by using generalized additive models (Klein et al., 2015; Rigby & Stasinopoulos, 2005; Umlauf & Kneib, 2018) or modern neural networks (Gasthaus et al., 2019; Rasp & Lerch, 2018). In related developments, semiparametric versions of quantile regression (Koenker, 2005) and transformation methods (Hothorn et al., 2014) can be leveraged for distributional regression.

Non-parametric approaches to distributional regression include kernel or nearest neighbour methods that depend on a suitable notion of distance on the covariate space. Then, the empirical distribution of the response for neighbouring covariates in the training set is used for distributional regression, with possible weighting in dependence on the distance to the covariate value of interest. Kernel smoothing methods and mixture approaches allow for absolutely continuous conditional distributions (Dunson et al., 2007; Hall et al., 1999; Li & Racine, 2008). Classification and regression trees partition the covariate space into leaves, and assign constant regression functions on each leaf (Breiman et al., 1984). Linear aggregation via bootstrap aggregation (bagging) or subsample aggregation (subagging) yields random forests (Breiman, 2001), which are increasingly being used to generate conditional predictive distributions, as proposed by Hothorn et al. (2004) and Meinshausen (2006).

Isotonicity is a natural constraint in estimation and prediction problems. Consider, for example, post-processing techniques in weather forecasting, where the covariates stem from the output of numerical weather prediction (NWP) models, and the response variable is the respective future weather quantity. Intuitively, if the NWP model output indicates a larger precipitation accumulation, the associated regression functions ought to be larger as well. Isotonic relationships of this type hold in a plethora of applied settings. In fact, standard linear regression analysis rests on the assumption of isotonicity, in the form of monotonicity in the values of the covariate(s), save for changes in sign.

Concerning non-parametric regression for a conditional functional, such as the mean or a quantile, there is a sizable literature on estimation under the constraint of isotonicity. The classical work of Ayer et al. (1955), Bartholomew (1959a, b), Brunk (1955), van Eeden (1958), Miles (1959) is summarized in Barlow et al. (1972), de Leeuw et al. (2009), Robertson et al. (1988). Subsequent approaches include Bayesian and non-Bayesian smoothing techniques (e.g. Dette et al., 2006; Mammen, 1991; Neelon & Dunson, 2004; Shively et al., 2009), and reviews are available in Groeneboom & Jongbloed (2014) and Guntuboyina & Sen (2018).

In distributional regression, it may not be immediately clear what is meant by isotonicity, and the literature typically considers one ordinal covariate only (e.g. Davidov & Iliopoulos, 2012; El Barmi & Mukerjee, 2005; Hogg, 1965; Rojo & El Barmi, 2003), with a notable exception being the work of Mösching & Dümbgen (2020b), whose considerations allow for a real-valued covariate. In the general case of a partially ordered covariate space, which we consider here, it is unclear whether semi- or non-parametric techniques might be capable of handling monotonicity contraints, and suitable notions of isotonicity remain to be developed.

To this end, we assume that the response Y is real-valued, and equip the covariate space $X$ with a partial order ≼. Our aim is to estimate the conditional distribution of Y given the covariate X, for short $L(Y|X)$⁠, on training data, in a way that respects the partial order, and we desire to use this estimate for prediction. Formally, a distributional regression technique generates a mapping from $x∈X$ to a probability measure $Fx$⁠, which serves to model the conditional distribution $L(Y|X=x)$⁠. This mapping is isotonic if $x≼x′$ implies $Fx≤stFx′$⁠, where $≤st$ denotes the usual stochastic order, that is, $G≤stH$ if G(y) ≥ H(y) for $y∈R$⁠, where we use the same symbols for the probability measures G, H and their associated conditional cumulative distribution functions (CDFs). Equivalently, $G≤stH$ holds if $G-1(α)≤H-1(α)$ for α ∈ (0, 1), where $G-1(α)=inf{y∈R:G(y)≥α}$ is the standard quantile function (Shaked & Shanthikumar, 2007).

Useful comparisons of predictive distributions are in terms of proper scoring rules, of which the most prominent and most relevant instance is the continuous ranked probability score (CRPS; Gneiting & Raftery, 2007; Matheson & Winkler, 1976). We show that there is a unique isotonic distributional regression that is optimal with respect to the CPRS (Theorem 1), and refer to it as the isotonic distributional regression (IDR). As it turns out, IDR is a universal solution, in that the estimate is optimal with respect to a broad class of proper scoring rules (Theorem 2). Classical special cases such as non-parametric isotonic quantile regression and probabilistic classifiers for threshold-defined binary events are nested by IDR. Simultaneously, IDR avoids pitfalls commonly associated with non-parametric distributional regression, such as suboptimal partitions of the covariate space and level crossing (Athey et al., 2019, p. 1167).

FIGURE 1

Simulation example for a sample of size n = 600 from the distribution in (1): (a) True conditional CDFs (smooth) and IDR estimates (step functions) for selected values of the covariate. (b) IDR estimated conditional distributions. The shaded bands correspond to probability mass 0.10 each, with the darkest shade marking the central interval. Vertical strips indicate the cross-sections corresponding to the values of the covariate in panel (a) [Colour figure can be viewed at https://dbpia.nl.go.kr/]

Open in new tabDownload slide

The remainder of the paper is organized as follows. The methodological core of the paper is in Section 2, where we prove existence, uniqueness and universality of the IDR solution, discuss computational issues and asymptotic consistency, and propose strategies for prediction. In Section 3, we turn to the critical issue of the choice of a partial order on the covariate space. Section 4 reports on a comparative simulation study that addresses both prediction and estimation, and Section 5 is devoted to a case study on probabilistic quantitative precipitation forecasts, with covariates provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) ensemble system. Precipitation accumulations feature unfavourable properties that challenge parametric approaches to distributional regression: The conditional distributions have a point mass at zero, and they are continuous and right skewed on the positive half-axis. In a comparison to state-of-the-art methods that have been developed specifically for the purpose, namely Bayesian model averaging (BMA; Sloughter et al., 2007), ensemble model output statistics (EMOS; Scheuerer, 2014) and heteroscedastic censored logistic regression (HCLR; Messner et al., 2014), the (out-of-sample) predictive performance of IDR is competitive, despite the method being generic, and being fully automatic once a partial order on the covariate space has been chosen.

We close the paper with a discussion in Section 6, where we argue that IDR provides a very widely applicable, competitive benchmark in probabilistic forecasting problems. The use of benchmark techniques has been called for across application domains (e.g. Basel Committee on Banking Supervision, 2016; Pappenberger et al., 2015; Rossi, 2013; Vogel et al., 2018), and suitable methods should be competitive in terms of predictive performance, while avoiding implementation decisions that may vary from user to user. IDR is well suited to this purpose, as it is entirely generic, does not involve any implementation decisions, other than the choice of the partial order, applies to all types of real-valued outcomes with discrete, continuous or mixed discrete-continuous distributions, and accommodates general types of covariate spaces.

2 Isotonic Distributional Regression

We proceed to introduce the IDR technique. To this end, we first review basic facts on proper scoring rules and notions of calibration. Then we define the IDR solution, prove existence, uniqueness and universality, and discuss its computation and asymptotic consistency. Thereafter, we turn from estimation to prediction and describe how IDR can be used in out-of-sample forecasting. Throughout, we identify a Borel probability measure on the real line $R$ with its cumulative distribution function (CDF), and we denote the extended real line by $R¯=[-∞,∞]$⁠.

2.1 Preliminaries

Introduced by Matheson & Winkler (1976), the CRPS has become popular across application areas and methodological communities, both for the purposes of evaluating predictive performance and as a loss function in estimation; see, for example, Gasthaus et al. (2019), Gneiting et al. (2005), Hersbach (2000), Hothorn et al. (2014), Pappenberger et al. (2015) and Rasp & Lerch (2018). The CRPS is reported in the same unit as the response variable, and it reduces to the absolute error, |x−y|, if F is the point or Dirac measure in $x∈R$⁠.

2.2 Existence, uniqueness and universality

A partial order relation ≼ on a set $X$ has the same properties as a total order, namely reflexivity, antisymmetry and transitivity, except that the elements need not be comparable, that is, there might be elements $x∈X$ and $x′∈X$ such that neither $x≼x′$ nor $x′≼x$ holds. A key example is the componentwise order on $Rd$⁠.

Let $I⊆R$ be an interval, and let S be a proper scoring rule with respect to a class $P$ of probability distributions on I that contains all distributions with finite support. Given training data in the form of a covariate vector $x=(x1,…,xn)∈Xn$ and response vector $y=(y1,…,yn)∈In$⁠, we may interpret any mapping from $x∈Xn$ to $Pn$ as a distributional regression function. Throughout, we equip $P$ with the usual stochastic order.

Definition 1

(S-based regression). An element $F̂=(F̂1,…,F̂n)∈Pn$ is an S-based isotonic regression of $y∈In$ on $x∈Xn$⁠, if it is a minimizer of the empirical loss

$ℓS(F)=1n∑i=1nS(Fi,yi)$

over all $F=(F1,⋯,Fn)$ in $P↑,xn$⁠.

In plain words, an S-based isotonic regression achieves the best fit in terms of the scoring rule S, subject to the conditional CDFs $F̂1,…,F̂n$ satisfying partial order constraints induced by the covariate values $x1,…,xn$⁠. The definition and the subsequent results can be extended to losses of the form $ℓS(F)=∑i=1nwiS(Fi,yi)$ with rational, strictly positive weights $w1,…,wn$⁠. The adaptations are straightforward and left to the reader.

Furthermore, the definition of an S-based isotonic regression as a minimizer of $ℓS$ continues to apply when $X$ is equipped with a pre- or quasiorder ≼ instead of a partial order. Preorders are not necessarily antisymmetric, and so there might be elements $x,x′$ such that $x≼x′$ and $x′≼x$ but $x′≠x$⁠. In this setting, we define x and $x′$ to be equivalent if $x≼x′$ and $x′≼x$⁠, and set $[x]≼p[x′]$ if representatives $u,u′$ of the equivalence classes $[x],[x′]$ satisfy $u≼u′$⁠. Then the binary relation $≼p$ defines a partial order on the set of equivalence classes, and the S-based isotonic regression with the new covariates and the partial order $≼p$ coincides with the original S-based isotonic regression.

In Supplementary Section S1 we prove the following result.

Theorem 1

(existence and uniqueness). There exists a unique CRPS-based isotonic regression  $F̂∈Pn$  of y on x.

At first sight, the specific choice of the CRPS as a loss function may seem arbitrary. However, the subsequent result, which we prove in Supplementary Section S1, reveals that IDR is simultaneously optimal with respect to broad classes of proper scoring rules that include all relevant choices in the extant literature. The popular logarithmic score allows for the comparison of absolutely continuous distributions with respect to a fixed dominating measure only and thus is not applicable here. Statements concerning calibration are with respect to the empirical distribution of the training data $(x1,y1),…,(xn,yn)$⁠.

Theorem 2

(universality) The IDR solution  $F̂$  of y on x is threshold calibrated and has the following properties.

• The IDR solution  $F̂$  is an S-based isotonic regression of y on x under any scoring rule of the form  

  $S(F,y)=∫(0,1)×RSα,θQ(F,y)dH(α,θ)$

  11

   or  

  $S(F,y)=∫R×(0,1)Sz,cP(F,y)dM(z,c),$

  12

   where  $Sα,θQ$  is the elementary quantile scoring function (6), $Sz,cP$  is the elementary probability scoring rule (7), and H and M are locally finite Borel measures on  $(0,1)×R$  and  $R×(0,1)$⁠, respectively.
• For every α ∈ (0, 1) it holds that  $F̂-1(α)=(F̂1-1(α),…,F̂n-1(α))$  is a minimizer of  

  $1n∑i=1nsα(θi,yi)$

  13

   over all  $θ=(θ1,…,θn)∈I↑,xn$⁠, under any function  $sα:I×I→R¯$  which is left-continuous in both arguments and such that  $S(F,y)=sα(F-1(α),y)$  is a proper scoring rule on  $P$⁠.
• For every threshold value z ∈ I, it is true that  $F̂(z)=(F̂1(z),…,F̂n(z))$  is a minimizer of  

  $1n∑i=1ns(ηi,1{yi≤z})$

  14

   over all ordered tuples  $η=(η1,…,ηn)∈[0,1]↓,xn$⁠, under any function  $s:[0,1]×{0,1}→R¯$  that is a proper scoring rule for binary events, which is left-continuous in its first argument, satisfies  $s(0,y)=limp→0s(p,y)$⁠, and is real-valued, except possibly s(0, 1) = −∞ or s(1, 0) = −∞.

The quantile weighted and threshold weighted versions of the CRPS studied by Gneiting & Ranjan (2011) arise from (11) and (12) with $H=G0⊗λ$ and $M=λ⊗G1$⁠, where λ denotes the Lebesgue measure, and $G0$ and $G1$ are σ-finite Borel measures on (0, 1) and $R$⁠, respectively. If $G0$ and $G1$ are Lebesgue measures, we recover the mixture representations (3) and (4) of the CRPS. By results of Ehm et al. (2016), if H is concentrated on ${α}×R$ and M is concentrated on {z } × (0,1), these representations cover essentially all proper scoring rules that depend on the predictive distribution F via $F-1(α)$ or F(z) only, yielding universal optimality in statements in parts (b) and (c) of Theorem 2.

In particular, as a special case of (13), the IDR solution is a minimizer of the quantile loss under the asymmetric piecewise linear or pinball function (5) that lies at the heart of quantile regression (Koenker, 2005). Consequently, as the mixture representation 2 of the CRPS may suggest, IDR nests classical non-parametric isotonic quantile regression as introduced and studied by Casady & Cryer (1976) and Robertson & Wright (1975). In other words, part (b) of Theorem 2 demonstrates that, if we (hypothetically) perform non-parametric isotonic quantile regression at every level α ∈ (0, 1) and piece the conditional quantile functions together, we recover the IDR solution. However, the IDR solution is readily computable (Section 2.3), without invoking approximations or truncations, unlike brute force approaches to simultaneous quantile regressions. Loss functions of the form (13) also include the interval score (Gneiting & Raftery, 2007, equation (43); Winkler, 1972), which constitutes the most used proper performance measure for interval forecasts.

In the special case of a binary response variable, we see from (c) and (14) that the IDR solution is an S-based isotonic regression under just any applicable proper scoring rule S. Furthermore, threshold calibration is the strongest possible notion of calibration in this setting (Gneiting & Ranjan, 2013, Theorem 2.11), so the IDR solution is universal in every regard. In the further special case of a total order on the covariate space, the IDR and pool adjacent violators (PAV) algorithm solutions coincide, and the statement in (c) is essentially equivalent to Theorem 1.12 of Barlow et al. (1972). In particular, the IDR or PAV solution yields both the non-parametric maximum likelihood estimate and the non-parametric least squares estimate under the constraint of isotonicity. The latter suggests a computational implementation via quadratic programming, to which we tend now.

2.3 Computational aspects

Clearly, a challenge in the computational implementation of IDR with general partial orders is that the number of variables in the quadratic programming problem (15) grows at a rate of $O(n)$⁠. As a remedy, we propose subsample aggregation, much in the spirit of random forests that rely on bootstrap aggregated (bagged) classification and regression trees (Breiman, 1996, 2001). It was observed early on that random forests generate conditional predictive distributions (Hothorn et al., 2004; Meinshausen, 2006), and recent applications include the statistical post-processing of ensemble weather forecasts (Schlosser et al., 2019; Taillardat et al., 2016; Taillardat et al., 2019). Bühlmann & Yu (2002) and Buja & Stützle (2006) argue forcefully that subsample aggregation (subagging) tends to be equally effective as bagging, but at considerably lower computational cost.

In view of the superlinear computational costs of IDR, smart uses of subsample aggregation yield major efficiency gains, taking into account that the estimation on different subsamples can be performed in parallel. Isotonicity is preserved under linear aggregation, and the aggregated conditional CDFs can be inverted to generate isotonic conditional quantile functions, with the further benefit of smoother estimates in continuous settings. A detailed investigation of optimal subsample aggregation for IDR is a topic for future research. For illustration, Figure 2 returns to the simulation example in Figure 1, but now with a much larger training sample of size n = 10,000 from the distribution in 1. Linear aggregation based on 100 subsamples (drawn without replacement) of size n = 1000 each is superior to the brute force approach on the full training set in terms of estimation accuracy. The computation on the full dataset for this simulation example takes 11.7 s for the naive implementation, but only 1.1 s for the sequential algorithm of Henzi et al. (2020). Subagging gives computation times of 9.9 and 2.5 s, respectively, or 1.8 and 0.5 s when parallelized over eight cores.¹

FIGURE 2

Simulation example for a sample of size n = 10,000 from the distribution in 1. The true conditional CDFs (smooth dashed graphs) are compared to IDR estimates (step functions) based on (a) the full training sample of size n = 10,000 and (b) linear aggregation of IDR estimates on 100 subsamples of size 1000 each [Colour figure can be viewed at https://dbpia.nl.go.kr/]

Open in new tabDownload slide

2.4 Consistency

We proceed to prove uniform consistency of the IDR estimator. While strong consistency of non-parametric isotonic quantile regression for single quantiles was proved decades ago (Casady & Cryer, 1976; Robertson & Wright, 1975), uniform consistency and rates of convergence for the IDR estimator have been established only recently, and exclusively in the case of a total order, see El Barmi & Mukerjee (2005, Theorem 1) and Mösching & Dümbgen (2020b, Theorem 3.3).

In Supplementary Section S2, we establish the following consistency theorem, which covers key examples of partial orders and is based on strictly weaker assumptions than the results of Mösching & Dümbgen (2020b). However, in contrast to their work, we do not provide rates of convergence. The choice $X=[0,1]d$ for the covariate space merely serves to simplify the presentation: As IDR is invariant under strictly isotonic transformations, any covariate vector $X=(X1,…,Xd)∈Rd$ can be transformed to have support in $[0,1]d$⁠, and the componentwise partial order can be replaced by any weaker preorder. A key assumption uses the concept of an antichain in a partially ordered set $(S,≼)$⁠, which is a subset $A⊆S$ that does not admit comparisons, in the sense that u ≼ v for u, v ∈ A implies u = v. As we discuss subsequently, results of Brightwell (1992) imply that the respective distributional condition is mild.

Theorem 3

(uniform consistency). Let  $X=[0,1]d$  be endowed with the componentwise partial order and the norm  $‖u‖=maxi=1,…,d|ui|.$  Let further  $(Xni,Yni)∈[0,1]d×R$⁠, $n∈N$⁠, $i=1,…,n$⁠, be a triangular array such that  $(Xn1,Yn1),…,(Xnn,Ynn)$  are independent and identically distributed random vectors for each  $n∈N$⁠, and let  $Sn={Xn1,…,Xnn}$⁠. Assume that

• (a)
  for all non-degenerate rectangles  $J⊆X$⁠, there exists a constant  $cJ>0$  such that  

  $#(Sn∩J)≥ncJ$

   with asymptotic probability one, that is, if  $An$  denotes the event that  $#(Sn∩J)≥ncJ$⁠, then  $P(An)→1$  as n → ∞;
• (b)
  for some γ ∈ (0,1),

  $max{#A:A⊂Snis antichain}≤nγ$

   with asymptotic probability one.

Assume further that the true conditional CDFs  $Fx(y)=P(Yni≤y|Xni=x)$  satisfy

• (c)

  $Fx(y)$  is decreasing in  $x$  for all  $y∈R$⁠;

• (d)
  for every  $η>0$⁠, there exists  $r>0$  such that  

  $sup{|Fx(y)-Fx′(y)|:x,x′∈[0,1]d,‖x-x′‖≤r,y∈R}<η.$

Then for every ε > 0 and δ > 0,

$limn→∞Psupx∈[δ,1-δ]d,y∈R|F̂x(y)-Fx(y)|≥ϵ=0.$

19

Assumption (a) requires that the covariates are sufficiently dense in $X$⁠, as is satisfied under strictly positive Lebesgue densities on $X$⁠. In order to derive rates of convergence, the size of the rectangles J in (a) would need to decrease with n, as in condition (A.2) of Mösching & Dümbgen (2020b); we leave this type of extension as a direction for future work. Assumption (c) is the basic model assumption of IDR, while assumption (d) requires uniform continuity of the conditional distributions, which is weaker than Hölder continuity in condition (A.1) of Mösching & Dümbgen (2020b).

Assumption (b), which is always satisfied in the case of a total order, calls for a more detailed discussion. In words, the maximal number of mutually incomparable elements needs to grow at a rate slower than $nγ$⁠. Evidently, the easier elements can be ordered, the smaller the maximal antichain. Consequently, Theorem 3 continues to hold under the empirical stochastic order and the empirical increasing convex order on the covariates introduced in Section 3.3, and indeed under any preorder that is weaker than the componentwise order. The key to understanding the distributional implications of (b) is Corollary 2 in Brightwell (1992), which states that for a sequence of independent random vectors from a uniform population on $[0,1]d$ the size of a maximal antichain grows at a rate of $n1-1/d$⁠; see also the remark following the proof of Theorem 3 in Supplementary Section S2.

As comparability under the componentwise order is preserved under monotonic transformations, any covariate vector $X∈Rd$ that can be obtained as a monotone transformation of a uniform random vector of arbitrary dimension guarantees (b). This includes, for example all Gaussian random vectors with non-negative correlation coefficients. In this light, assumption (b) is rather weak, and well in line with the intuition that for multivariate isotonic (distributional) regression to work well, there ought be at least minor positive dependence between the covariates. In the context of our case study in Section 5, high positive correlations between the covariates are the rule, as exemplified by Table 3 in Raftery et al. (2005).

2.5 Prediction

As noted, the IDR solution $F̂=(F̂1,…,F̂n)∈P↑,xn$ is defined at the covariate values $x1,…,xn∈X$ only. Generally, if a (not necessarily optimal) distributional regression $F=(F1,…,Fn)∈P↑,xn$ is available, a key task in practice is to make a prediction at a new covariate value $x∈X$ where $x∉{x1,…,xn}$⁠. We denote the respective predictive CDF by F.

The difference between the bounds (if any) in (20) might be a useful measure of estimation uncertainty and could be explored as a promising avenue towards the quantification of ambiguity and generation of second-order probabilities (Ellsberg, 1961; Seo, 2009). In the context of ensemble weather forecasts, the assessment of ambiguity has been pioneered by Allen & Eckel (2012). Interesting links arise when the envelope in (20) is interpreted in the spirit of randomized predictive systems and conformal estimates as studied by Vovk et al. (2019); compare, for example, their Figure 5 with our Figure 4b below.

3 Partial Orders

The choice of a sufficiently informative partial order on the covariate space is critical to any successful application of IDR. In the extreme case of distinct, totally ordered covariate values $x1,…,xn∈X$ and a perfect monotonic relationship to the response values $y1,…,yn$⁠, the IDR distribution associated with $xi$ is simply the point measure in $yi$⁠, for i = 1, …, n. The same happens in the other extreme, when there are no order relations at all. Hence, the partial order serves to regularize the IDR solution.

Thus far, we have simply assumed that the covariate space $X$ is equipped with a partial order ≼, without specifying how the order might be defined. If $X⊆Rd$⁠, the usual componentwise order will be suitable in many applications, and we investigate it in Section 3.1. For covariates that are ordinal and admit a ranking in terms of importance, a lexicographic order may be suitable.

If groups of covariates are exchangeable, as in our case study on quantitative precipitation forecasts, other types of order relations need to be considered. In Sections 3.2 and 3.3, we study relations that are tailored to this setting, namely, the empirical stochastic order and empirical increasing convex order. Proofs are deferred to Supplementary Section S3.

3.1 Componentwise order

This order becomes weaker as the dimension d of the covariate space increases: If $x~=(x1,…,xd,xd+1)$ and $x~′=(x1′,…,xd′,xd+1′)$ then $x≼x′$ is a necessary condition for $x~≼x~′$⁠. The following result is an immediate consequence of this observation and the structure of the optimization problem in Definition 1.

Proposition 1

Let  $x=(x1,…,xn)$  and  $x*=(x1*,…,xn*)$  have components  $xi=(xi1,…,xid)∈Rd$ and $xi*=(xi1,…,xid,xi,d+1)∈Rd+1$ for $i=1,…,n$⁠, and let S be a proper scoring rule. Then if  $Rd$  and  $Rd+1$  are equipped with the componentwise partial order, and  $F̂$  and  $F̂*$  denote S-based isotonic regressions of y on x and  $x*$⁠, respectively, it is true that  

$ℓS(F̂*)≤ℓS(F̂).$

In simple words, under the componentwise partial order, the inclusion of further covariates can only improve the in-sample fit. This behaviour resembles linear regression, where the addition of covariates can only improve the (unadjusted) R-square.

3.2 Empirical stochastic order

We now define a relation that is based on stochastic dominance and invariant under permutation.

Definition 2

Let  $x=(x1,…,xd)$  and  $x′=(x1′,…,xd′)$  denote elements of  $Rd$⁠. Then x is smaller than or equal to  $x′$  in empirical stochastic order, for short  $x≼stx′$⁠, if the empirical distribution of  $x1,…,xd$  is smaller than the empirical distribution of  $x1′,…,xd′$  in the usual stochastic order.

This relation is tailored to groups of exchangeable, real-valued covariates. The following result summarizes its properties and compares to the componentwise order, which we denote by ≼.

Proposition 2

Let  $x=(x1,…,xd)$  and  $x′=(x1′,…,xd′)$  denote elements of  $Rd$  with order statistics  $x(1)≤⋯≤x(d)$  and  $x(1)′≤⋯≤x(d)′$⁠.

• (a)

  The relation  $x≼stx′$  is equivalent to  $x(i)≤x(i)′$ for i=1, …, d.

• (b)

  If  $x≼x′$ then $x≼stx′$⁠.

• (c)

  If  $x≼stx′$  and x and  $x′$  are comparable in the componentwise partial order, then  $x≼x′$⁠.

• (d)

  If  $x≼stx′$  and  $x′≼stx$  then x and  $x′$  are permutations of each other. Consequently, the relation  $≼st$  defines a partial order on  $R↑d$⁠.

In a nutshell, the empirical stochastic order is equivalent to the componentwise order on the sorted elements, and this relation is weaker than the componentwise order. However, unlike the componentwise order, the empirical stochastic order does not degenerate as further covariates are added. To the contrary, empirical distributions of larger numbers of exchangeable variables become more informative and more easily comparable.

3.3 Empirical increasing convex order

In applications, the empirical stochastic order might be too strong, in the sense that it does not generate sufficiently informative constraints. In this light, we now define a weaker partial order on $R↑d$⁠, which also is based on a partial order for probability measures. Specifically, let X and $X′$ be random variables with CDFs F and $F′$⁠. Then F is smaller than $F′$ in increasing convex order if $E(ϕ(X))≤E(ϕ(X′))$ for all increasing convex functions ϕ such that the expectations exist (Shaked & Shanthikumar, 2007, Section 4.A.1).

Definition 3

Let  $x=(x1,…,xd)$  and  $x′=(x1′,…,xd′)$  denote elements of  $Rd$⁠. Then x is smaller than or equal to  $x′$  in empirical increasing convex order, for short  $x≼icxx′$⁠, if the empirical distribution of  $x1,…,xd$  is smaller than the empirical distribution of  $x1′,…,xd′$  in increasing convex order.

This notion provides another meaningful relation for groups of exchangeable covariates. The following result summarizes its properties and relates it to the empirical stochastic order.

Proposition 3

Let  $x=(x1,…,xd)$  and  $x′=(x1′,…,xd′)$  denote elements of  $Rd$  with order statistics  $x(1)≤⋯≤x(d)$  and  $x(1)′≤⋯≤x(d)′$⁠.

• (a)
  The relation  $x≼icxx′$  is equivalent to  

  $∑i=jdx(i)≤∑i=jdx(i)′forj=1,…,d.$
• (b)

  If  $x≼stx′$  then  $x≼icxx′$⁠.

• (c)
  If  $x≼icxx′$  then  

  $1d∑i=1dxi+d-12(d+1)g(x)≤1d∑i=1dxi′+d-12(d+1)g(x′),$

   where g is the Gini mean difference,  

  $g(x)=1d(d-1)∑i,j=1d|xi-xj|.$

  22

• If  $x≼icxx′$  and  $x′≼icxx$ then x and $x′$  are permutations of each other. Consequently, the relation  $≼icx$  defines a partial order on  $R↑d$⁠.

Figure 3 illustrates the various types of relations for points in the positive quadrant of $R2$⁠. As reflected by the nested character of the regions, the componentwise order is stronger than the empirical stochastic order, which in turn is stronger than the empirical increasing convex order. The latter is equivalent to weak majorization as studied by Marshall et al. (2011). In the special case of vectors with non-negative entries, their Corollary C.5 implies that $x∈Rd$ is dominated by $x′∈Rd$ in empirical increasing convex order if, and only if, it lies in the convex hull of the points of the form $(ξ1xπ(1)′,…,ξdxπ(d)′)$⁠, where π is a permutation and $ξi∈{0,1}$ for i = 1, …, d.

FIGURE 3

Regions of smaller, greater and incomparable elements in the positive quadrant of $R2$⁠, as compared to the point (1, 3), for the (left) componentwise, (middle) empirical stochastic and (right) empirical increasing convex order. Coloured areas below (above) of (1, 3) correspond to smaller (greater) elements, while blank areas contain elements incomparable to (1, 3) in the given partial order [Colour figure can be viewed at https://dbpia.nl.go.kr/]

Open in new tabDownload slide

4 Simulation Study

Since we view IDR primarily as a tool for prediction, we compare it to other distributional regression methods in terms of predictive performance on continuous and discrete, univariate simulation examples, as measured by the CRPS. However, as noted below and formalized in Supplementary Section S4, the CRPS links asymptotically to $L2$ estimation error, so under large validation samples, prediction and estimation are assessed simultaneously. A detailed comparative study on mixed discrete-continuous data with a multivariate covariate vector is given in the case study in the next section.

Implementation details for the various methods are given in Table S1 in Supplementary Section S5. Here we only note that QRF uses the grf package (Tibshirani et al., 2020) with a splitting rule that is tailored to quantiles (Athey et al., 2019). We see that, unlike IDR, its competitors rely on manual intervention and tuning. For example, QRFs perform poorly under the default value of 5 for the tuning parameter min.node.size, which we have raised to 40. Further improvement may arise when tuning parameters, such as honesty fraction and node size, are judiciously adjusted to the specific scenario and training sample size at hand. In contrast, IDR is entirely free of implementation decisions, except for the subagging variant, $IDRsbg$⁠, where we average predictions based on estimates on 100 subsamples of size n/2 each.

Table 1 shows the mean CRPS for the different methods and simulation scenarios. Scenario (23) is the same as in the introduction and illustrated in Figure 1. It has a smooth covariate–response relationship, and NP, SQR, and even the misspecified TRAM technique, which are tailored to this type of setting, outperform QRF and IDR. However, the assumption of continuity in the response is crucial, as the results under the discontinuous scenario (24) demonstrate, where IDR and $IDRsbg$ perform best. In the non-isotonic scenario (25), IDR and $IDRsbg$ retain acceptable performance, even though the key assumption is violated. Not surprisingly, SQR faces challenges in the Poisson scenario (26), where the conditional quantile functions are piecewise constant, and IDR is outperformed only by TRAM. Throughout, the simplistic subagging variant of IDR has slightly lower mean CRPS than the default variant that is estimated on the full training set, and it would be interesting to explore the relation to the super-efficiency phenomenon described by Banerjee et al. (2019).

These results lend support to our belief that IDR can serve as a universal benchmark in probabilistic forecasting and distributional regression problems. For sufficiently large training samples, IDR offers competitive performance under any type of type of linearly ordered outcome, without reliance on tuning parameters or other implementation choices, except when subsampling is employed.

5 Case Study: Probabilistic Quantitative Precipitation Forecasts

The past decades have witnessed tremendous progress in the science and practice of weather prediction (Bauer et al., 2015). Arguably, the most radical innovation consists in the operational implementation of ensemble systems and an accompanying culture change from point forecasts to distributional forecasts (Leutbecher & Palmer, 2008). An ensemble system comprises multiple runs of numerical weather prediction (NWP) models, where the runs or members differ from each other in initial conditions and numerical-physical representations of atmospheric processes.

Ideally, one would like to interpret an ensemble forecast as a random sample from the conditional distribution of future states of the atmosphere. However, this is rarely advisable in practice, as ensemble forecasts are subject to biases and dispersion errors, thereby calling for some form of statistical post-processing (Gneiting & Raftery, 2005; Vannitsem et al., 2018). This is typically done by fitting a distributional regression model, with the weather variable of interest being the response variable, and the members of the forecast ensemble constituting the covariates, and applying this model to future NWP output, to obtain conditional predictive distributions for future weather quantities. State of the art techniques include Bayesian Model Averaging (BMA; Raftery et al., 2005; Sloughter et al., 2007), Ensemble Model Output Statistics (EMOS; Gneiting et al., 2005; Scheuerer, 2014), and Heteroscedastic Censored Logistic Regression (HCLR; Messner et al., 2014).

In this case study, we apply IDR to the statistical post-processing of ensemble forecasts of accumulated precipitation, a variable that is notoriously difficult to handle, due to its mixed discrete-continuous character, which requires both a point mass at zero and a right skewed continuous component on the positive half-axis. As competitors to IDR, we implement the BMA technique of Sloughter et al. (2007), the EMOS method of Scheuerer (2014), and HCLR (Messner et al., 2014), which are widely used parametric approaches that have been developed specifically for the purposes of probabilistic quantitative precipitation forecasting. In contrast, IDR is a generic technique and fully automatic, once the partial order on the covariate space has been specified.

5.1 Data

The data in our case study comprise forecasts and observations of 24 h accumulated precipitation from 6 January 2007 to 1 January 2017 at meteorological stations on airports in London, Brussels, Zurich and Frankfurt. As detailed in Table 2, data availability differs, and we remove days with missing entries station by station, so that all types of forecasts for a given station are trained and evaluated on the same data. Both forecasts and observations refer to the 24 h period from 6:00 UTC to 6:00 UTC on the following day. The observations are in the unit of millimetre and constitute the response variable in distributional regression. They are typically, but not always, reported in integer multiples of a millimetre.

Statistical post-processing is both a calibration and a downscaling problem: Forecasts and observations are at different spatial scales, whence the unprocessed forecasts are subject to representativeness error (Wilks, 2019, Chapter 8.9). Indeed, if we interpret the predictive distribution from the raw ensemble (27) as the empirical distribution of all 52 members—a customary approach, which we adopt hereinafter—there is a strong bias in probability of precipitation forecasts: Days with exactly zero precipitation are predicted much less often at the NWP model grid box scale than they occur at the point scale of the observations.

5.2 BMA, EMOS and HCLR

Before describing our IDR implementation, we review its leading competitors, namely, state of the art parametric distributional regression approaches that have been developed specifically for accumulated precipitation.