Breiman's "Two Cultures" Revisited and Reconciled

\noindentKeywords: Integrated statistical learning theory, Exploratory machine learning, Uncertainty prediction machine, ML-powered modern applied statistics, Information theory. \vskip.8em {.2ex} { \setcounter{tocdepth}{2} } \vskip.5em \linespread{1.28} {1.4ex}

Introduction

A Clash of Two Cultures

By the dawn of the 21st century, the statistics community was clearly divided into two distinct camps--namely, the (parametric) statistics and the (algorithmic) machine learning. Leo breiman01 presents a vivid description of this cultural polarization in his influential essay on “Statistical Modeling: The Two Cultures.” But what has triggered this split? \vskip.05em ML$>$STAT: breiman01 argued that algorithmic models are far more flexible, scalable, and accurate for complex big data problems. The traditional statistical methods, by contrast, are based on a priori assumed parametric models that are mainly suitable for small datasets.

\vskip.05em STAT$>$ML: This algorithmic-supremacy standpoint had been furiously debated by several eminent statisticians in the discussion of the original paper. They found Breiman's claim problematic on several counts. First, despite convincing prediction results, the absence of an explicit data-generating model (a probabilistic base) can make algorithmic methods less useful for scientific investigation; see the commentaries by cox01 and parzen01. efron2020pred articulated this admirably: “Abandoning mathematical models comes close to abandoning the historic scientific goal of understanding nature.” Second, the black-box nature of algorithmic models makes them inscrutable, opaque, and not easily interpretable. In fact, regarding machine learning models, Breiman himself agreed that “their mechanism for producing a prediction is difficult to understand...So on interpretability, they rate an F.” \vskip.05em Where do we stand now? To a large extent, the cultural division between these two statistical learning frameworks has been growing steadily over the past decades, especially in the post-deep-learning era. What options are we left with? 1) Let's get rid of black-boxes and use statistical parametric models, even though that means sacrificing accuracy; 2) Let's stick to black-boxes because they are accurate, even though that means sacrificing interpretability and robustness. Perhaps unsurprisingly, neither of these two extreme positions is tenable. To make real progress, we are forced to consider the question: how do these two different cultures fit together? Can we build a framework of data analysis that combines the best ideas from both sides and offer some way to make two into one? In this paper, we show that such an integrated framework of statistical learning is entirely buildable.

The message of this paper is very simple: the algorithmic models should not be feared because of their complexity, parametric models should not be shunned because of their simplicity, we can connect them both to create a more powerful learning scheme. With the help of a synthetic example, the next section explains what we mean by `powerful' and why we need such a technology.

The Butterfly Prediction Problem

We are given $(x_i,y_i)$, $i=1,\ldots,N$ as shown in the Fig (ref), which we refer to the butterfly data. The goal is to build a prediction model for the outcome variable $Y$ given $X=x$.

The conventional regression-based prediction methods proceed as follows: \vskip.25em (1) Estimation: We have a vast number of choices available. On the one hand, we have the simplest parametric model--the linear regression--and on the other, a large collection of sophisticated machine learning (ML) methods like k-nearest neighbors (knn), random forest, support-vector machine (svm), gradient boosting machine (gbm), neural net, etc. Irrespective of whether we choose a parametric or algorithmic method, the estimated-regression function takes the form, more or less, of a `flat' zero line, as shown in Fig. (ref). \vskip.2em (2) Inference: The obvious next question is, how accurate is the estimate? The blue dotted line in Fig. (ref) displays the 90% bootstrap confidence band, which is narrow at the center (around $0$) and gets wider as we move outward--correctly reflecting the uncertainty of the mean function. Conclusion: the small bootstrap error (uncertainty) indicates that the estimated regression line is impressively accurate! \vskip.25em {\bf But, how useful is this analysis?} Even though we have an “accurate” estimate of the underlying regression prediction, it does not tell us anything useful. In fact, it paints a grossly misleading picture of the reality (in the sense that $x$ is not useful for predicting $y$, which is clearly untrue). So, what went wrong? Do we need a bigger, more complex model? Better tuning of the hyperparameters? More computing power? More data? Unfortunately, `more data + more computing power' does not necessarily guarantee a better algorithm. The problem lies elsewhere.

\vskip.1em {\bf How to avoid statistical sin in the practice of regression?} To avoid “elegantly solving the wrong problem” we must first ask: what to estimate, before tackling the question of how to estimate (linear or nonlinear; parametric or algorithmic; Bayes or frequentist). Almost all prediction methods start with the following `$\mu+\ep$' model\,\footnote{Open your favorite machine learning textbook, and it won’t be long before you encounter `$\mu+\ep$' model.}, which approximates $y$ as a non-linear function of $x=(x_1,\ldots,x_p)$: \beq y = \mu(x)\,+\,\ep, \eeq where the random errors $\epsilon_i$'s are i.i.d with expected value zero, and finite variance; often assumed to be Gaussian with fixed variance $\ep \sim \cN(0,\si^2)$. The conventional statistical learning estimates the high-dimensional mean (regression) function $\mu(x)$, which can then be used for prediction. Now, it is not hard to see that no matter how complex the ML method we choose to approximate the function $\mu(x)$, it will still be fruitless (and a turbo-waste of computation\footnote{Algorithms are not bulldozer, they are tools for understanding.}) since the conditional mean contains no useful information for the butterfly example! But then, how should we perform prediction for the butterfly data? What should we report? These are basic questions, yet often overlooked in the rush to apply fancy ML methods.

\vskip.1em {\bf Uncertainty prediction machine}. The main object of interest is the random variable $Y|X=x$, which will be denoted by $Y_{|x}$. Can we really predict the value of the response variable $Y$ for a given $X=x$? No. It impossible to foresee what exact value the random variable $Y_{|x}$ will take\footnote{After all, Heisenberg (1927) taught us, we can never know enough about the present to exactly predict the future event, except for predicting probabilities of different possible outcomes--the indeterminacy principle.}. At best, we can predict the probabilities and its distribution over all possible values of the response variable--the conditional distribution $f_{Y|X=x}(y)$, which encodes a complete description of $Y_{|x}$. The problem of statistical learning is to design an `uncertainty prediction machine' to infer $f_{Y|X=x}(y)$ from the observed data $(x_i,y_i)$, $i=1,\ldots,N$. Classical ML methods generally predict some type of summary statistic (e.g., location or scale) of the conditional distribution. However, the butterfly example taught us that this limited information could be dangerously misleading when dealing with messy heterogeneous datasets. One of the motivating questions behind our study is: how can we convert “first generation” machine learning methods into uncertainty distribution prediction machine? Our interest lies in extracting the underlying statistical “law” (governing equation) that is driving the data.

Uncertainty Prediction Machine: What to Expect?

Before delving into the technical details, it may be worthwhile if we start by taking a `crude look at the whole' to clearly convey what practical results can be hoped from our statistical learning technology. Given the data $(x_1,y_1),\ldots,(x_N,y_N)$ with $x_i \in \cR^p$ and $y_i \in \cR$, our primary focus is the inference of $Y_{|x}$, where $Y_{|x}$ denotes the conditional random variable $Y|X=x$ with distribution function $F_{Y|X=x}(y)$ and density $f_{Y|X=x}(y)$. We will be using the butterfly data as a stylized example to describe the different stages of our analysis.

\vskip.24em {\bf Level 0}. ML curve-fitting. We start with a user-selected machine learning method ML$_0$ to estimate (possibly high-dimensional) regression function $\widehat{\mu}(x)$ along with its uncertainty band. The result for butterfly data is already displayed in Fig. (ref). Frankly speaking, this is where the practice of traditional machine learning ends, and ours begins.

\vskip.24em {\bf Level 1}. Heterogeneity diagnostic. A natural question is whether the initial ML$_0$ captured the essential relationship between $y$ and the predictor variables. This can be checked by searching for (unmodeled) patterns in the residuals. Recall that for the butterfly data $\hat \mu(x)\approx 0$ for all $x$; thus residuals are actually the original response values $y_1,\ldots,y_n$.

The critical task boils down to checking whether the residuals are homogeneous, i.e., $f_{Y|X=x}(y)$ $=f_Y(y)$ for all $x$. This would ensure that ML$_0$ was able to capture all the important predictive information. However, in most practical real-life problems (including butterfly data) the residuals will be heterogeneous. And the main challenge is to identify the components of $f_{Y|X=x}(y)$ that are significantly changing with $x$. We like to address this by constructing a simple diagnostic plot that can help users to visually examine and detect the prominent sources of heterogeneity in the residuals. Fig. (ref)(a) is one such graphical exploratory plot drawn for the butterfly data, which reveals scale (second-order effect) and tail (fourth-order kurtosis effect) of $f_{Y|X=x}(y)$ are changing significantly with $x$. Accordingly, this tool helps us to discover the hidden patterns that, while invisible to ML$_0$, still persist and should be modeled for sharpening the prediction quality.

\vskip.24em

{\bf Level 2}. Pivot density. We now turn our attention to the problem of predicting conditional distribution at a particular $x$, say at $x_0=2$ for the butterfly example. There are some obvious candidates: $\wtf_Y(y)$ the empirical (marginal) distribution of $y$; or its nonparametric smooth (say, kernel-smoothed) version $\fhat_Y(y)$; or some kind of parametric model motivated from domain-knowledge or statistical convenience, e.g., $\cN(\hat \mu(x_0),\hat \si_y^2)$, where $\hat \si_y$ denotes the standard error of $y$. We call these initial starting “guesses” the pivot density; they are nothing but a crude first-approximation of the true unknown conditional density estimate. For the butterfly data, we choose $\cN(0,\hat \si_y^2=2.55^2)$ to be the pivot, as shown in Fig. (ref)(b).

A natural question to ask is whether the presumed pivot model is approximately correct. In most practical cases, the answer will be no, they are not adequate. Thus we have to take some corrective measures to rectify the deficiency of the starting pivot. But before doing so, it is crucial to devise some technique that can reveal the deficiency of the current model. The only problem is that we don't have enough $y$-samples at $X=x_0$ to perform traditional goodness-of-fit tests. We have to come up with something new that is computable and interpretable.

\vskip.24em {\bf Level 3}. Uncertainty quantification for the pivot. The basic question is whether the starting pivot density $f_0(y)$ is a good model for $Y_{|x_0}$. We answer this question by estimating the following contrast density function for $u=F_0(y)$ \beq d(u;F_0,F_{Y|X=x}):=\,d_x^0(u)\,=\dfrac{f_{Y|X=x} ( F_0^{-1}(u))}{f_0( F_0^{-1}(u))}, 0<u<1\eeq where $F_0^{-1}(u)$ is the quantile function for the pivot $f_0(y)$. Note that $d_x^0(u)$ captures the deviation of $f_0(y)$ from the true unknown conditional density $f_{Y|X=x_0}(y)$, and by virtue of doing so it helps us to quantify and characterize the uncertainty of the pivot density. By inspecting the shape of $d_{x_0}^0(u), 0<u<1$ (how it deviates from uniformity), practitioners can quickly infer whether $f_0(y)$ needs any repair or not.

Fig. (ref)(c) displays the estimated contrast function for the butterfly data. Looking at the shape, it becomes clear that the initial Gaussian pivot completely missed the bimodality of $f_{Y|X=2}(y)$. The contrast function $d_x^0(u)$ plays a central role in our theory. In Section (ref) we provide an estimation strategy. The fascinating aspect of the algorithm is that it utilizes the power of the initial machine learning algorithm ML$_0$ to get the covariate-adaptive $\whd_x^0$.

\vskip.24em {\bf Level 4.} Density modulation. When the contrast function indicates the inadequacy of the starting pivot $f_0(y)$ as a model for the conditional distribution $f_{Y|X=x}(y)$, what to do next? How to modify (or `boost') $f_0(y)$ to make it data-consistent? This can be solved in a remarkably simple way via density modulation: $\hf_{Y|X=x}=f_0 \times \whd_x^0$, where $d_{x}^0$ modulates (modifies) the shape of starting $f_0(y)$ to produce the conditional density. The detailed theory of why and how this works will be discussed in Section 2. Fig. (ref)(d) displays the estimated predictive density for the butterfly data, which we can write down analytically as \beq \hf_{Y|X=x}(y),=\,\dfrac{1}{\hat \si} \phi\big(\frac{y}{\hat \si}\big) \times \big\{1 -0.20\, T_2(y;F_0) - \,0.47\,T_4(y;F_0)\big\}, where $\hat \si=2.55$\eeq which is product of two components: the Gaussian pivot and the contrast density $d_{x_0}^0(F_0(y))$. The basis functions $\{T_j(y;F_0)\}$ in equation (ref) are specially-designed rank-based orthonormal polynomials of the pivot. They injects robustness and non-linearity into the data modeling process--but more on this later. We will also show that this new class of $d$-modulated conditional density models can generate an extraordinarily rich spectrum of shapes.

One notable aspect of (ref) is that we have `identified' a parametric model (instead of blindly assuming, which Breiman sharply criticized) by algorithmic (fully nonparametric) means. In Section (ref), we elucidate the process by which we boost and convert the “$0$th order” machine learning method (ML$_0$) into a finitely parameterizable uncertainty prediction machine. This allows our learning philosophy to integrate the fundamental character of both algorithmic data-driven modeling and parametric statistical modeling. Nevertheless, the predicted distribution (ref) compactly represents all possible outcome values along with their respective probabilities, which is essential for decision-making under uncertainty.

\vskip.24em {\bf Level 5}. Predictive inference: What are some of the most likely values of the response variable $y$, given that we have observed $X=x_0$? For a fixed $\al$, ideally, we would like to find the smallest volume region that covers $100(1-\al)\%$ area of the conditional distribution, aka the highest density prediction region. The orange highlighted part in Fig. (ref)(d) displays 68% highest-density prediction region for the butterfly data, which consists of two disjoint intervals [-3.33,-1.40] $\cup$ [1.20,3.28]. The volume of this set (the sum of the lengths of the intervals) can be used to quantify the uncertainty of the prediction. Additionally, we can extract all the conditional quantiles from the estimated $\whF_{Y|X=x}(y)$; see panel (e). Our quantile regression curves are, by design, non-crossing and nonlinear. Curious readers may contrast this with the traditional quantile regression estimator, as shown in Fig. (ref) of Appendix A.4.1.

\vskip.24em {\bf Level 6}. Verification: Data modeling typically consists of four stages: estimation (of the initial pivot), exploration (discovering the deficiency of the pivot), rectification (repairing the pivot), and verification (checking whether the corrected pivot emulates the observed reality). The last step remains to be done, which requires evidence that the estimated uncertainty prediction model provides a satisfactory description of the data. Scientific data analysis, by definition, needs to be testable (falsifiable), and there's no two ways about it.

For our butterfly example, we compute $U_i=\whF_{Y|X=x_i}(y_i)$ for $i=1,\ldots,N$ using our trained model on a 15% hold-out validation dataset. It is not difficult to see that the distribution of $U_i$'s will be `close' to Uniform[0,1], if the estimated conditional density can explain the essential patterns/variations in the data. The quantiles of the empirical distribution of $U_i$'s are plotted against the respective quantiles of the uniform distribution at the bottom right of Fig. (ref). The tight cluster of points around the 45-degree reference line strongly indicates that the model is successful in capturing the underlying data generating process, and can be trusted to make predictions. Additionally, in Section (ref), we provide a formal goodness-of-fit test for uniformity, which, when applied for this example, yields a p-value of $0.82$.

Goals and Contributions

Figure (ref) shows a schematic description of the modular architecture of our data analysis pipeline, which aims to “assemble” different statistical learning cultures--algorithmic, parametric, exploratory, and robust modeling--into one coherent whole.

$\bullet$ Unified Interface: We like to construct a generic “interface” that can convert any algorithmic regression procedure into an uncertainty prediction machine. We are less interested in solutions that depend on the specific choice of an ML method; developing `tricks' on a case-by-case basis is too clumsy. This raises the question: how to design such a simple and unified interface? Section (ref) presents some general theory that enables such a design.

\vskip.1em $\bullet$ Algorithmic Parameterization: The next crucial aspect of our method lies in distilling an explicit and simple statistical model from a black-box technique. Section (ref) lays out such a mechanism to allow algorithmic synthesis of probabilistic models with interpretable statistical parameters. This opens up the possibility of building closed-form, interpretable statistical models that are as expressive and scalable as machine learning models. \vskip.1em

$\bullet$ Exploratory Machine Learning: How can we check whether a trained machine learning method is consistent with the observed data? If not, can we discover its blind spots using some exploratory graphical tools? Finally, can we revise the starting algorithmic method when required? All these basic questions are often left unanswered by traditional machine learning procedures. In the words of Richard hamming97art, “AI has traditionally stuck to the what is done and seldom considered the how it is done.” The `how' part is important to build trusted AI/ML systems. Section (ref) discusses some new tools to make the `modeling process' more transparent and interpretable.\footnote{ A recently-passed legislation of New York City Council says that “As we advance into the 21st century, we must ensure our government is not black boxed.”} It consists of the following chain of well-connected steps: \vskip.1em Data + Pivot $\rightarrow$ Exploration $\rightarrow$ Discovery $\rightarrow$ Rectification $\rightarrow$ Verification $\rightarrow$ Prediction. \vskip.05em This layout of data analysis embodies the viewpoint that empirical (statistical) science is a `self-correcting' process.

\vskip.15em $\bullet$ `Robust + Efficient' Machine Learning: Current ML systems are extremely vulnerable to small perturbations of the training data, which makes them ineffective for critical applications. One of the unsolved dilemmas is to construct learning algorithms that are simultaneously robust against noise and leaks very little, if any, efficiency under the standard (unperturbed) situation. Section (ref) develops a new robust learning theory to tackle this important issue. \vskip.2em

$\bullet$ Fundamentals of Applied Statistics: The new mechanics of data modeling that is presented here is not a “one-trick pony,” it has deep connections with the basic fundamental statistics. In section (ref), we demonstrate how one can systematically derive and generalize (to large high-dimensional problems in a scalable and flexible manner) many traditional and modern statistical methods as a special case of our framework. This could be especially useful for the practice and pedagogy of ML-powered modern applied statistics. Ultimately, interlinking Statistics and Machine Learning will enrich and revitalize both the communities, by creating excitement to work across the boundary.

Integrated Statistical Learning Theory

In this section, we introduce the basic theoretical underpinnings, principles, and methods of the integrated statistical-learning framework.

Entropy, Uncertainty, and Prediction

The marginal distribution $f_Y(y)$ embodies uncertainty of $Y$. Now, once we observe $X=x$, the uncertainty is updated to $f_{Y|X=x}(y)$. We say that some information is received through the observation $X=x$, if the conditional distribution $f_{Y|X=x}(y)$ differs from the initial $f_Y(y)$. Therefore, one can measure the additional information gain on $Y$ after observing $X=x$ by examining the change in the probability distributions $f_{Y}(y)$ and $f_{Y|X=x}(y)$. Information theory provides an elegant framework for quantifying this “change.”

{\bf Proof}: First substitute $F_Y(y)=u$ in (ref), and then apply the definition (ref) to deduce: \beq {\rm KL}(d_x;\mathbbm{1}_{[0,1]}) = \int_0^1 d_x(u) \log d_x(u) \dd u.\eeq

For each $x_i$ in our dataset, we will get a different $d_{x_i}$ and a different KL-divergence value (ref). Can we measure the overall predictive power of $X$ by taking average of all these different KL-divergence values? To satisfactorily answer that, we have to introduce the concept of entropy.

Predictability-index can be defined as the difference between the entropy of $Y$ and $Y|X$, which measures how much of our uncertainty about $Y$ decreased after observing $X$. The following theorem establishes a beautiful connection between predictability-index and our conditional information density function.

The proof is given in Appendix A.3. The important point here is that predictive information can be expressed solely in terms of $d_x$, which distills the discrepancy (the “gap”) between the marginal $f_Y$ and the conditional $f_{Y|x}$. Accordingly, conditional information density (CID) $d_x(u)$ acts as a natural starting point for developing a statistical theory of predictive modeling.

$d$-Modulated Conditional Density Representation

A few key points about this conditional density formula:

$\bullet$ This new representation of conditional density is universally valid and does not impose any distributional assumptions on the outcome (e.g., dispersion, symmetry, skewness, heavy-tailedness, etc.). Most importantly, it decouples the conditional density into two components: the marginal $f_Y$ (the `pivot' density) and the conditional information density $d_x$ (which encapsulates all the essential knowledge for predicting $Y$ from $X=x$).

$\bullet$ We call the density perturbation formula of Eq. (ref) the $d$-modulation of marginal $Y$. The function $d_x$ filters out the directions by which the initial $f_Y(y)$ departs most significantly from the true $f_{Y|X=x}(y)$. For that reason, we call it, alternatively, the contrast density function.

We now provide another interpretation of $d_x(u)$ by introducing the notion of conditional probability integral transform (shortly, cPIT).

{\bf Proof}: $D_x(u)=\Pr(F_Y(Y_{|x}) \le u)=F_{Y|x}( F_Y^{-1}(u))$. Taking derivative we get the density: \[ D'_x(u) \,\equiv\, d_x(u)\,=\, \dfrac{f_{Y|x}(F_Y^{-1}(u)) }{f_Y(F_Y^{-1}(u))},~0<u<1.\]

A Robust Learning Theory

It is evident from the decomposition (ref) that we need to estimate $d_x(F_Y(y))$ to obtain the conditional density of $Y$ given $X=(X_1,\ldots,X_p)$. In this section, we describe a special technique for estimating $d_x$ that enjoys the following desirable properties:

$\bullet$ Robust as well as efficient: It has been already noticed that (i) existing statistical learning methods are notoriously sensitive to noise, which makes them an easy target for adversarial attacks; and at the same time, (ii) naively constructed robust algorithms perform poorly in the ideal setting with no contamination. The real question comes down to this: how to robustify algorithms so that it can be resistant to noise/outliers while also making sure it is sufficiently accurate (efficient) when applied to relatively clean data.\footnote{This `robustness-efficiency' dilemma (also known as `rank-and-size' duality) is not new. It is at least $60$ years old unsettled issue of Statistics; see box1953non and tukey1960.} Here we will present a new class of statistical learning methods to balance robustness and efficiency.

$\bullet$ Scalable to massive high-dimensional data: An important practical goal is to have a method that works for large-$(n,p)$ problems. We achieve this goal by utilizing the power of general machine learning tools.

$\bullet$ Interpretable: It would be nice to have a tractable (and preferably smooth) model for the contrast density $d_x(u)$ with a few interpretable parameters that can capture increasingly complex shapes of the conditional distributions.

Its a challenging task to reliably learn the “true” relationship between the input features $X$ and output $Y$ from big noisy datasets. To ensure that we succeed at it, our method has to be `doubly-robust,' and thus tackle possible contamination in both response and covariates. We will approach the estimation of $d_x$ in two stages: first, we will make it robust with respect to the outcome variable $Y,$ and then we will discuss defence strategies to guard against corrupted input $X$.

{\bf $Y$-Robustness}. Substitute $u=F_Y(y)$ in the definition (ref), to see that $d_x(F_Y(y))$ allows us to represent the ratio $f_{Y|X=x}(y)/f_Y(y)$ as a function of the rank-transform (i.e., probability integral transform) $F_Y(Y)$. This is the first major step towards ensuring the $Y$-robustness of our method. As a result, one can expand $d_x(F_Y(y))$ in the orthonormal basis of $F_Y(y)$. One such orthonormal system is the LP-family of rank-polynomials Deep17LPMode,deep16MetaLP,D20copula,DeepLPKsample19, which we denote $\{T_j(y;F_Y)\}$ to emphasize that they are polynomials of $F_Y(Y)$ not $Y$, hence extremely robust. As the true $F_Y$ is unknown, we will instead use the empirical LP-bases (eLP) $\{T_j(y;\wtF_Y)\}$ for our data analysis. Appendix A.2 describes its construction and properties.

Summing up the above discussion, we have the following LP-series representation: \beq d_x(F_Y(y)) = 1\,+\,\sum_{j} \LP_{j|x} T_j(y;F_Y). \eeq The problem of estimating the function $d_x$ now boils down to estimating the orthogonal coefficients $\LP_{j|x}$. To get an expression for the conditional LP-coefficients, first recall that $T_j$'s are mean zero random variables and are orthonormal with respect to measure $F_Y$, i.e., \beq \int T_j(y;F_Y) T_k(y;F_Y) f_Y(y) \dd y = \delta_{jk}, \eeq where $\delta$ is the Kronecker delta function. Consequently, we multiply both sides of the equation (ref) with $T_j(y;F_Y)\cdot f_Y(y)$ and then integrate to produce the following identity: \beq \LP_{j|x} = \int_{-\infty}^\infty d_x(F_Y(y))\, T_j(y;F_Y)\,f_Y(y) \dd y, for $j=1,2,\ldots$ \eeq This immediately yields the following important theorem.

The proof is trivial; just substitute $d_x( F_Y(y))$ in equation (ref) by $f_{Y|X=x}(y)/f_Y(y)$. Theorem (ref) allows us to interpret the conditional coefficients $\LP_{j|x}$ (which are surely function of $x$) as a specialized regression function, where we regress $T_j(y;F_Y)$ on $X$. This is extremely handy, because one can now use any general machine learning method to estimate these coefficients!

\vskip.34em {\bf $X$-Robustness}. Nevertheless, the computational theory (ref) is still not fully satisfactory, as corrupted (or adversarially perturbed) $X$ can substantially influence the coefficients. To provide further defense, we need one last piece of idea, which is a little-known but fundamental fact about `regression via rank-transform.' Also see Appendix A.4.2 for something surprising.

{\bf Proof}: We start with a basic quantile mechanics fact: for any random variable $X$ we have \[Q_X(F_X(X))~=~X, ~~\text{with probability 1}.\] where $Q_X(u)=\inf_x\{F_X(x) \ge u\}$, $0<u<1$ denotes the quantile function of $X$. Accordingly, any function of $X$, say $h(X)$, can be equivalently rewritten as a function of the rank-transform $F_X(X)$, since $h(X)$ is equal to $h\circ Q_X(F_X(X))$ with probability $1$.

In plain language, this result proves the “sufficiency” of rank-transform (retaining all the information) for estimating conditional mean function. Therefore, we can estimate $\LP_{j|x}$ (ref) by regressing $T_j(Y;\wtF_Y)$ on the LP-bases of $X$. This will allow doubly-robust estimation of $d_x$ without sacrificing an iota of efficiency in the process. Next, we will go over some examples, to demonstrate how this theory actually translates into practice (algorithm). \vskip.35em

{\bf Example 1.} The Butterfly data. Fig. (ref) displays the scatter plots of $T_j(y,\wtF_Y)$ against $u=\wtF_X(x)$ for $j=1,\ldots,4$. We refer these graphs as LP-scatter plots. The red curves denote knn-estimated regression curves. In the R computing language, this is simply knn($T_j(y;\wtF_Y)$\, $\sim$ $\,T_X$), $(j=1,\ldots,4)$ where $T_k(x;\wtF_X)$ is the $k$th col of LP-transformed feature matrix $T_X$. Of course, instead of knn method, one can choose any machine learning routine; the next example uses gradient boosting algorithm (gbm). To get the estimated conditional LP-coefficients at $x_0=2$ (indicated by the blue dots) we record the predicted values at $u_0=\wtF_X(x_0)$: $\widehat\LP_{2|x_0}$ equals to $-.020$, $\widehat\LP_{4|x_0}$ equals to $-0.47$, and the rest are practically zero. Substituting them into (ref) yields the estimated contrast function $\whd_{x_0}$, as was shown in Fig. (ref)(c).

{\bf Example 2.} BUPA liver disorders data. This is a popular dataset (available in the UCI ML repository) used by many researchers, including breiman01. It consists of measurements of gamma-glutamyl transpeptidase (GGT) and alanine-aminotransferase (ALT) extracted from $N=345$ male individuals' blood samples. High levels of ALT and GGT indicate a damaged liver. Fig. (ref) (a) shows the scatter plot of $y$=log(ALT) versus $x$=log(GGT). We seek to predict the distribution of log-ALT levels given $x=3.5$; see the red triangle in panel (a) of Fig. (ref).

Our data analysis scheme can be simply summarized as follows:

1) We start by modeling the conditional mean function using gradient-boosting algorithm. The rugged blue line denotes the estimated regression curve $\widehat{\mu}(x)$.

2) To estimate the conditional density at $x_0=3.5$, we choose $\cN(\widehat{\mu}(x_0), \si_y^2)$ as our pivot model, where $\widehat{\mu}(x_0)=3.35$ and $\si_y=0.5$.

3) Following theorems 4 and 5, we estimate the conditional LP-Fourier coefficients $\LP_{j|x_0}$ and the associated contrast density function: $\whd_{x_0}(F_0(y))~ \approx ~1 - 0.72\, T_2(y;F_0) + 0.10\, T_3(y;F_0).$ This is shown in the middle panel of Fig. (ref).

4) Interestingly, we can interpret the coefficients of $\whd_{x_0}$: the large value of $\widehat \LP_{2|x_0}$ indicates that the selected pivot needs scale-correction (“second-order” correction) to get close to the true conditional density at $x_0$. In addition, the negative sign implies the variability has to be reduced, which is vividly apparent from Fig. (ref)(c).

5) Similarly, the presence of third-order (positive) LP-coefficient indicates the symmetric normal pivot needs to be tilted to the right. The positive skewness of $f_{Y|x_0}(y)$ is probably justified due to the presence of few large $y$-values\footnote{Classical robust methods treat them as “outliers” and devise mechanisms for removing them from analysis. This is incorrect. A careful examination will reveal that these large-$y$ values (indicated by squares) have some pattern; they are not randomly scattered (like outliers). The real question of what is an outlier and what is a “real” pattern, finally boils down to artfully balancing robustness with efficiency.} around $x_0=3.5$, as indicated by squares.

6) Finally, the shape of $\whd_x$ provides quick exploratory guidance for choosing appropriate families of parametric conditional distributions. For the BUPA liver-disorders data, one could select skew-normal or asymmetric logistic distribution FriedmanCD2020 as a model for $Y_{|x}$.

The Basic Algorithm

We are now in a position to provide the blueprint of the algorithmic interface in Fig. (ref) that converts a user-selected machine learning method into an uncertainty prediction machine in a robust and interpretable manner. The theory in the previous section, makes it absolutely straightforward and extremely easy to implement, which essentially involves three steps. The following algorithm presents a flowchart of the core computational steps in pseudo-code, so that one can implement it using any programming language or ML-platform. We have found that the H2O open-source is extremely efficient. But users can choose any ML environment (e.g., TensorFlow, Azure, AutoGloun, etc.) to operationalize the following workflow.

\hrule height .65pt \vskip.25em {\bf Input} Data $(X_i,y_i)$ for $i=1,\ldots, N$; user-selected machine learning method (ML); values of $m_y$ and $m_x$. And the target test-point $x_0 \in \cR^p$.

Step 1. $T_Y \leftarrow$ LP.basis$(y,m_y)$; $T_X \leftarrow$ LP.basis$(X,m_x)$

for($j=1,\ldots,m_y$)

\{

Step 2. Fitted.ML$_j$ $\leftarrow$ ML\big($T_Y[\,, \,j] \sim T_X$\big)

Step 3. $ \LP_{j|x_0} \leftarrow$ predict(Fitted.ML$_j$, $x_0$)

\}

{\bf Return} $\whd_{x_0} \,\leftarrow \,1 + \sum_j \LP_{j|x_0} T_Y[\,, \,j]$

\hrule height .65pt \vskip1em

The LP-transformation is implemented by the function LP.basis; $T_X$ is the stacked matrix $[T_{X_1}\mid T_{X_2} \mid \cdots \mid T_{X_p}]$, where the $j$th column of the matrix $T_{X_\ell}$ is simply $T_j(x;\wtF_\ell)$. Step 1 makes the whole procedure robust, non-linear, and invariant under monotone transformations. Step 2 performs the training of the ML method and step 3 predicts the fitted value at the desired $X=x_0$. Based on our experience, $m_x=4$ and $m_y=6$ work reasonably well for most problems, including all the datasets analyzed in this paper. Finally, one can generate the uncertainty distribution $f_{Y|x_0}$ by applying the $d$-modulation formula of equation (ref).

$d$-Kernel Smoothing and Importance Weighting

In many applications, we are often interested in estimating quantities like $\Ex[\Psi(Y)|X=x]$ for general $\Psi(\cdot)$ function. This can be rewritten as \[\Ex[\Psi(Y)|X=x]\,=\int_y \Psi(y) f_{Y|X=x}(y)\dd y\,=\int_y \Psi(y) f_Y(y) \,d_x(F_Y(y)) \dd y,\] by applying the formula (ref). This leads to the following theorem, which expresses the (local) conditional mean as a weighted average of the global samples $\{y_1,\ldots,y_N\}$.

Exploratory Learning Machine

Our goal here is to develop some graphical diagnostic methods to make a transition from a black-box answer machine to an exploratory learning machine that can drill deeper and generate previously unanticipated, interesting questions to facilitate data-driven scientific discovery and understanding.\footnote{A `better-question' is often better than a `better-answer.'} This can be viewed as a pre-deployment “health check-up” for statistical learning methods that comes with a set of nonparametric exploratory tests.

Heterogeneity Component Analysis

Heterogeneity component analysis (HCA) is a technique for identifying dynamic components of $Y_{|x}$--the aspects of the conditional distribution $f_{Y|X=x}(y)$ (such as location, scale, skewness, tail, etc.) that vary with covariates. We illustrate the main idea using the butterfly example. The procedure starts with LP-scatter plots, as was shown in Fig. (ref). The first scatter plot captures how the location of $f_{Y|X=x}(y)$ (conditional mean) changes with $x$; the second one captures the change of scale (the nature of heteroscedasticity); the third one provides information on the changing skewness and so on. Thus, an obvious approach for detecting the heterogeneity components would be to:

Step 1. Compute a goodness-of-fit statistic for the $j$th LP-scatter plot by its coefficient of multiple correlation: \[R^2_j\, =\, \text{Proportion of variance in $T_j(y;F_Y)$ that is explained by $\{T_1(x;F_X), \ldots, T_m(x;F_X)\}$}\] Step 2. Assess the significance of the $R^2_j$ by constructing F-statistic for the $j$th LP-scatter \[ F_j = \dfrac{R^2_j}{1-R^2_j}\times\dfrac{N-m-1}{m}\] whose null distribution\footnote{To allow non-normal errors, use the fact that $m F_{m,N-m-1}$ approaches to $\chi^2_m$ (in the sense of convergence in distribution) as $N \rightarrow \infty$ and compute the p-value based on this large-sample chi-square approximation.} is $F_{m,N-m-1}$. Fig. (ref)(a) displays the significant (whose pvalue $<0.05$) heterogeneity components for the butterfly data. \vskip.25em {\bf The “Flatland” Hypothesis}. For real data analysis, one can apply the above method on residuals to determine which shape parameters are varying with $x$, beyond first-order mean. The “flatness” of residuals implies the errors are homogeneous with respect to the covariates.

{\bf Example 3.} Bone mineral density data hastie1999bone. It contains measurements on the relative change in spinal bone density (over one year period) of $N=485$ North American adolescents, as a function of age. This dataset was previously analyzed by stanfybook. Fig. (ref)(a) shows the data along with the estimated regression smoother. We like to understand whether there is any “excess” heterogeneity beyond the 1st-order regression modeling. To check the homogeneity of the residuals with respect to the covariate $x$, we apply our HCA procedure on the residuals $y_i-\widehat y_i$. The result is shown in the top row of Fig. (ref), which says that the regression model has not taken into account the heteroscedasticity of the data (furthermore, there is some evidence of varying skewness in the conditional distribution). These unmodelled (dynamic) heterogeneity components should be modeled to improve the prediction. The heterogeneity is partly due to the presence of both male and female youngsters' in the dataset. \vskip.15em

\vskip.35em {\bf Example 4.} Online news popularity data fernandes2015. The objective of this study is to predict the popularity of online news. The data\footnote{The article in row \# $31,038$ seems to be an outlier (in fact it is an erroneous data: check its 4th and 5th feature values--rate of unique and non-stop words--which can't be more than 1, by definition). Since our method is doubly-robust we don't worry that much. But standard machine learning methods might get unduly influenced by this one data point (article). The practice of looking at the data is always a good habit.} consist of $N = 39,644$ articles, and for each article we have $y$ (log$_{10}$ of number of shares in social networks, a measure of popularity) and $p=59$ extracted features (such as number of words in the title, number of links, number of images, etc.) $x$. This dataset was recently analyzed by FriedmanCD2020.

We start by fitting gradient boosting regression to the data. At this point, it is a natural question to ask whether the gbm captured most of the important information in the data. If not, what is missing? To spot the unmodelled dynamic components, we start from the residuals, as shown in Fig. (ref)(c). The HCA diagnostic plot of panel (d) is constructed based on the Lasso-selected features, which clearly indicates the presence of a strong heteroscedasticity. However, it is also true that the conditional distribution will be asymmetric and long-tailed, but our analysis shows that they are not swiftly changing with covariates; they are more or less static in nature, not dynamically evolving shape-parameters. This further corroborates the findings of FriedmanCD2020.

Pivot Uncertainty Modeling

The parameters $\LP_{j|x}$ of the contrast density $d_x$ can be viewed as “generalized coordinates” that describe the shape of $f_{Y|X=x}(y)$ relative to some reference distribution (pivot) $f_0(y)$.

In fact, these shape-coefficients $\LP_{j|x}$ capture the whole dynamics of how the conditional distribution $f_{Y|X=x}(y)$ evolves with $x$.\footnote{Indeed, one can even construct an $m$-dimensional phase-space (analogous to statistical mechanics), consisting of $N$ points $\{\LP_{1|x_i},\ldots,\LP_{m|x_i}\}_{i=1}^N$ to graphically represent the evolving shape of $f_{Y|X=x_i}(y)$. Each point in LP-phase-space corresponds to a particular shape for a particular value of $x=(x_{1},\ldots,x_{p})$.} Naturally, a large value of \beq \operatorname{qPivot}(x)\,=\, \sum_{j}\big| \LP_{j|x} \big|^2 \eeq indicates higher uncertainty (misfit) of the pivot $f_0$ for modeling the conditional $f_{Y|X=x}$. The `q' in $\operatorname{qPivot}$ stands for quantification of lack-of-fit. The above formula can also be interpreted as measuring the deviation of $d_x(u)$ from uniformity (follows from Parseval's identity): \beq\int_0^1 (d_x(u)-1)^2 \dd u = \sum_{j} \big| \LP_{j|x} \big|^2 \eeq

{\bf Example 5.} Movie box-office revenue data filmdata2012. The goal of this study is to build a forecasting model for film revenues. The data contain $N=4031$ pairs of observations $(x_i,y_i)$ on the 1990s film data, where $x_i$ is the log of opening box-office revenues and $y_i$ is the log of box office revenues after the first week. We are interested in predicting the distribution of the outcome $y$ for $x=12$ and $14$--marked with the blue and green dotted line in Fig. (ref) (a). Our starting pivot $f_0 =\cN(14,\si_y=3.25)$ is shown in the panel (b). To check whether the conditional density at $x$ follows this presumed parametric law, we estimate the respective contrast densities; see the middle panel of Fig. (ref). The bimodality of $\whd_{x=12}$ indicates that the Gaussian pivot failed to capture the presence of films with two kinds of box-office earnings. This is also reflected in the pattern of the scatter plot, which is divided into two branches. On the other hand, the inverted `U' shape (quadratic) of $\whd_{x=14}$ implies that the variability of the Gaussian pivot needs to be adjusted ($2$nd-order correction) to go from $f_0(y) \mapsto f_{Y|X=14}(y)$. The final refurbished conditional density models at $x=12$ and $14$ are shown in the bottom row of Fig. (ref). They are estimated using the $d$-modulation technique, described in Sec. (ref).

Goodness-of-fit Diagnostics

The flexibility of our approach (as depicted in Fig. (ref)) enables us to choose any machine learning method to design the uncertainty prediction machine (UPM). Nonetheless, a user may want to know which ML-powered distribution prediction algorithm(s) better explains the pattern in the data. To answer that here we introduce a graphical exploratory procedure to evaluate the overall goodness-of-fit of the method. Deploying statistical models without knowing how well it actually fits the data, is a statistical sin. We illustrate our idea using the butterfly data. \vskip.25em

{\bf ML-driven UPM system}. To construct an UPM for the butterfly data, we start with three possible base learners: k-nearest neighbors (knn), random forest (RF), and gradient-boosting algorithm (GBM). We train each ML-driven UPM model using the generic procedure described in section (ref). To check whether these models are congruent with the observed data, we proceed as follows:

Step 1. Compute the generalized quantile-residuals for all three models on a hold-out validation dataset (separate from the training samples) of size $n_v$ \beq U_{\ell i} = \whF_{Y|X=x_i}^\ell(y_i), for $i=1,\ldots,n_v$; \,and $\ell=1,2,3$. \eeq Note that if the model accurately captures the pattern in the validation data, then we expect the distribution of the $U_{\ell}$'s to be close to uniform. \vskip.2em

Step 2. Display the histograms and QQ-plot of $\{U_{\ell 1}, \ldots, U_{\ell n_v}\}$. This is done in Fig. (ref) for the butterfly data, which shows the knn model satisfactorily captures the observed reality.

\vskip.2em

Step 3. To quantify the model-fit we propose the following test statistic for the $\ell$-th model:

\beq \operatorname{qDIV}_\ell = \sum_{j=1}^m \Big | n_{v}^{-1} \sum_{i=1}^{n_v} \Leg_j \big(U_{\ell i} \big) \Big|^2, \eeq where $\Leg_j$ denotes the $j$th orthonormalized shifted Legendre polynomial on unit interval. The $\operatorname{qDIV}$ statistic measures the distance between Uniform$[0,1]$ and the sample distribution of $U_{\ell}$. We have used $m=6$ in all our examples. One can also compute pvalues using the $\chi^2_m$ null distribution of $\operatorname{qDIV}$. For more details see Deep18DMT and DeepLPKsample19.

\vskip.2em

Step 4. It is also interesting to note that one can directly compute these generalized quantile-residuals from the contrast distribution function, since (by substituting $F_0(t)=u)$: \[ F_{Y|X=x}(y)\,= \int _{-\infty}^y f_0(t) \,d_x(F_0(t))\dd t\,=\int_0^{F_0(y)} d_x(u) \dd u \,= \,D_x(F_0(y)).~~ \vspace{-.15em} \]

For the butterfly data $\operatorname{qDIV}$ values (with p-values in the parentheses) for knn, random forest, and gbm respectively are $0.0466 \,(0.82)$, $1.12 \,(1.4\times10^{-10})$, and $0.386\, (0.0026)$. These statistic values can be used as a general measure of “lack-of-fit” for comparing (ranking or even tuning hyperparameters) the performances of different ML-engine based UPMs.

Sharpening Black-box Generators

Imagine we are given some black-box generative model that can simulate $y$-samples for a given $x=(x_1,\ldots,x_p)$: \[\texttt{MODEL}(x)\, ~\mapsto \,~\{\wty^x_1,\ldots, \wty^x_s \}.~~~\] Without having any information on the particular MODEL\footnote{This could be any arbitrarily complex computational program, e.g., deep learning method with millions of parameters, an ensemble of $100$ machine learning models, or the method proposed in FriedmanCtree2019.}, our goal is to: (i) check whether the samples are truly generated from the underlying conditional distribution $f_{Y|X=x}(y)$; (ii) if not, explain `why' by displaying the associated $\whd_x$; and finally, (iii) sharpen the given imperfect samples to make them compatible with the underlying stochastic data generating mechanism.

\vskip.3em {\bf Butterfly data}. To illustrate the method, let us consider that we are given all the $y_i$'s for which the $x_i$ value is between $0$ and $4$. The histogram of these samples is shown in the left panel of Fig. (ref). The blue curve shows the true $f_{Y|X=2}(y)$. Starting from these weak conditional samples, our goal is to refine it so that they comply with the underlying law at $x=2$. \vskip.15em Step 1. Choose the empirical distribution $\widecheck{f}_{y|x}$ of $\{\wty^x_1,\ldots, \wty^x_s \}$ be the pivot; see Fig. (ref)(a). \vskip.15em Step 2. Randomly generate one sample $\wty_{\mbox{*}}$ from $\widecheck{f}_{y|x}$ and $U$ from Uniform$[0, 1]$. \vskip.15em Step 3. Let $\widecheck{F}_{y|x}$ denotes the cdf of $\widecheck{f}_{y|x}$. Then, set $\widehat{y}^x=\wty_{\mbox{*}}$, if \[d_x\big( \widecheck{F}_{y|x} ( \wty_{\mbox{*}}) \big)~>~ U \max_u d_x(u). \] Otherwise discard it and go back to the previous step. Here $d_x(u)$ means $d_x(u;\widecheck{F}_{y|x},F_{y|x})$. \vskip.15em Step 4. Repeat steps 2-3 until we have $s$-sample $\{\widehat{y}^x_1,\ldots, \widehat{y}^x_s\}$, which we call $d$-sharp samples, since $d_x$ acts as a tool for sharpening the initial weak conditional samples; see Fig. (ref)(b).

Some Connections and Unified Interpretation

The theoretical constructs of our integrated learning framework are not based on an isolated “trick,” but are fundamental to statistics. Our formalism and principles of data modeling unify many traditional and modern statistical methods. In this section, we will report findings to highlight this fact.

K-sample Comparison Problems

Given random samples from $k$ distributions $F_i(y),$ $i=1,\ldots,k$ we seek to test the following null hypothesis: \beq H_0: F_1(y) = \cdots = F_k(y) for all y, \eeq with the alternative hypothesis that at least two distributions are different. Before attacking this general problem of equality of distributions, we start with a simpler problem of testing the equality of means. The Kruskal–Wallis statistic provides a nonparametric test for it.

Definition and notation. We have $k$ mutually independent random samples with sizes $n_1,\ldots n_k$ with the combined sample size $N$; $R_i$ denotes the rank of $Y_i$ in the pooled samples; $G_\ell$ is the set of indices for the $\ell$-th group. The Kruskal–Wallis test statistic is defined as: \beq {\rm KW} = \dfrac{12}{N(N+1)} \sum_{\ell=1}^k n_\ell \Big\{ \frac{1}{n_\ell} \sum_{i \in G_\ell} R_i - \frac{N+1}{2} \Big\}^2. \eeq One may go one step further and ask how to test the equality of scale parameters (variability over different groups), which can be tested using Mood statistic, defined as \beq {\rm Mood} = \dfrac{180}{N(N+1)(N^2-4)} \sum_{\ell=1}^k n_\ell\Big\{ \frac{1}{n_\ell} \sum_{i \in G_\ell} \big(R_i - \frac{N+1}{2}\big)^2 - \frac{N^2-1}{12} \Big\}^2 \eeq For testing high-order differences like skewness (not to mention kurtosis, which will take almost two full lines to write), the formula of the test statistic gets even more complicated.\footnote{ If you are a machine learning practitioner and concerned that these are too abstract and complicated, hang in there--don't give up. We will soon see they can be represented in a beautiful and compact manner using our modern notation, which can be implemented in one line, without the burden of memorizing any formula.} {\beq {\rm SKEW}=\frac{7}{N(N+1)(N^2-4)(N^2-9)}\sum_{\ell=1}^k n_\ell\Bigg\{ \frac{20}{n_\ell} \sum_{i \in G_\ell} \big(R_i - \frac{N+1}{2}\big)^3 - \dfrac{3N^2-7}{n_\ell} \sum_{i \in G_\ell} \big( R_i - \frac{N+1}{2}\big)\Bigg\}^2 \eeq \vskip.025em } A practitioner at this point might wonder how to make sense of these monstrous formulae? Is any intuitive or systematic derivation possible? The answer is yes. We just have to look at the problem from a different perspective. A modern introduction to the $k$-sample comparison problem is given below. \vskip.15em Step 1. View it as a $(X,Y)$ problem: Define $X$ to be the group-index variable taking values $\{1,\ldots, k\}$. We now have $\{(x_i,y_i)\}_{i=1}^N$ data, which can be displayed as a scatter plot. \vskip.35em Step 2. Reformulate $H_0$: The $k$-sample null-hypothesis of equality of distributions can now be viewed as a test of the homogeneity problem (see Section (ref)) $f_{Y|X=\ell}(y)=f_Y(y)$ for $\ell=1,\ldots,k$. Now since \beq f_{Y|X=\ell}(y)\,=\,f_Y(y) \Big\{ 1\,+\sum_j \LP_{j|\ell} T_j(y;F_Y) \Big\},\eeq the original $k$-sample hypothesis (ref) can alternatively be rewritten as \[H'_0:~ \LP_{j|\ell}=0, ~~\text{for all}~ j.~~~\] Hence, intuitively, it makes sense to apply the HCA diagnostic test of section (ref) on the $(x,y)$ data. So what happens, if we dare to take the leap of faith? We get almost all of the known results of $k$-sample modeling in one-shot.\footnote{After all, “the man who perpetually hesitates accomplishes nothing; It is the man who dares who wins.”} To justify the statement, we start with the following theorem.

The proof is in Appendix A.3. The beauty of this result is that it seamlessly generalizes to high-order cases\footnote{The agreement of our general theory with the $k$-sample results may appear to be a `miraculous coincidence' but in reality, this actually reveals the fundamental nature of it.} (proofs are not difficult and are left to the reader): \beas MOOD&=&\sum_{\ell=1}^k n_\ell \,\big| \tLP_{2|\ell} \big|^2\,=\,N R_2^2; (change in scale)\\ SKEW&=& \sum_{\ell=1}^k n_\ell \,\big| \tLP_{3|\ell} \big|^2\,=\,N R_3^2; (change in skewness)\\ KURT&=& \sum_{\ell=1}^k n_\ell \,\big| \tLP_{4|\ell} \big|^2\,=\,N R_4^2, (change in tail) \eeas and so on. The practical benefit of these results is that they provide a systematic and comprehensive $k$-sample analysis program that is radically simple to implement, since LP-multiple correlations $R_j^2$ can be computed in a single line R-code: \texttt{summary}(\texttt{lm}($T_j(y) \sim T_X$)\$\texttt{r.squared}. \vskip.2em {\bf Example 6.} LDL cholesterol of Quail hettmansperger2010book. This is a randomized experiment investigating LDL (low-density lipoprotein) cholesterol in quails; $N=39$ quails were randomly assigned to $k=4$ diets for a specific period of time, where each diet mixed with a different drug compound. The boxplot is shown in Fig. (ref). The scientific question is: whether or not different drug compounds change the LDL cholesterol levels. The right panel of Fig. (ref) shows the re-scaled HCA-plot, where we multiply the $j$th LP-multiple-correlations $R_j^2$ by the sample size $N$. It shows that the distribution of LDL cholesterol levels is changing in scale and location.

Probabilistic Index Model

Do baseball players gain weight as they get older? Can a larger dose of an antipsychotic drug reduce the severity of depression (after adjusting for the confounders)? Do working-class Belgian families spend less of their income on food, proportionately, as they come to make more money? All of these problems are concerned with modeling the comparison probability regression: \beq {\rm CPR}(x) = \Pr(Y \le Y_{|x}), \eeq as a function of $x$, instead of as a usual conditional mean regression. {\rm CPR}$(x)$ encodes the probability that a randomly selected subject from the population with $X=x$ has a higher response value than a randomly selected subject from the full data. This concept was pioneered by Olivier Thas and his collaborators under the name of probabilistic index model; see thas2012PIM. To understand how (ref) is related to our framework, we start with the following integral representation: \beq \Pr(Y \le Y_{|x}) = \int F_Y(y) \dd F_{Y|X=x}(y).\eeq This leads to the following important result.

The first equality follows from (ref), and the second equality follows from supplementary equation (ref). The representation (ref) justifies the name comparison probability regression, which is simply a standardized 1st-order LP-regression function. As a consequence of Theorem (ref), the whole procedure can be implemented in one line of R-code: lm$(T_1(y)\sim T_X)$.

\vskip.3em {\bf Example 7}. Baseball data matloff2017MLbook. Fig. (ref) shows the age verses weight scatter plot of $N=1015$ major-league baseball players. We would like to investigate whether baseball players gain weight as they age. In other words, our task is to check whether $\Pr(Y \le Y_{|x})$ increases as a function of age or not. We proceed as follows:

Step 1. Compute the LP-basis functions: $T_1(y;\wtF_Y)$ for the response and $T_X=[T_1(x;\wtF_X),$ $\ldots,T_m(x;\wtF_X)]$ for the predictor variable. Here we have picked $m=4$.

Step 2. (penalized) LP-regression: Perform the multiple linear regression $T_1(y;\wtF_Y)$ on $T_X$. Return the Bayesian information criterion (BIC) selected model. If no variable is selected then it supports the null hypothesis $H_0: \Pr(Y \le Y_{|x})=0.5$. For the baseball data we get: \beq \tLP_{1|x}=\Ex[T_1(y;\wtF_Y)|X=x]=0.17\,T_1(x;\wtF_X). \eeq There is no intercept as $\Ex[T_1(y;\wtF_Y)]=0$, by construction. Although athletes strive to keep physically fit, eq. (ref) seems to suggest that the baseball players do gain some weight (in a slow linear rate) over time. Fig. (ref)(b) makes this somewhat clear to understand how.

Shape Predictors: Generalized Feature Selection

When dealing with a large number of covariates, it is often important to understand which are the most predictive for the response variable $Y$. Currently, the most popular feature selection method is $L_1$-penalized lasso regression stanfybook, which minimizes the residual sum of squares $\| y-X\beta\|^2_2$ subject to $\|\be\|_1 \le \lambda$. This will be denoted algorithmically as lasso$(y\sim X)$. It is important to note that, when we perform lasso-regression of $y$ on $X$, we only get `first-order' informative features--the variables that influence the conditional mean $\Ex[Y|X=x]$. Naturally, the question is: can we go beyond mean, and find the most relevant attributes to describe the shape of the conditional distribution $f_{Y|X=x}(y)$? Can we better utilize the lasso-technique to find those “shape predictors”? The following is an attempt to construct such a generalized feature selection program, which is extremely easy to implement. The algorithmic logic proceeds as follows:

Step 1. Let's start by recalling our basic model: $f_{Y|X=x}(y)=f_Y(y)\cdot d_x(F_Y(y))$. The covariates $x=(x_1,\ldots,x_p)$ can influence the distribution of the response $Y$ only through $d_x \circ F_Y(y)=1 + \sum_{j}\LP_{j|x}T_j(y;F_Y)$, which is governed by the following system of equations: \beas \LP_{1|x}&=&\Ex[T_1(Y;F_Y)|X=x]\\ \LP_{2|x}&=&\Ex[T_2(Y;F_Y)|X=x]\\ \vdots & & \vdots \eeas The first equation captures the dynamics of conditional mean as a function of $x$, while the second one describes how the variability (scale) changes with $x$, and so on.

Step 2. Thus, it is quite obvious that one can extract the $j$th order features, denoted by $$\bO_j\,=\,\big\{\ell: \widehat \be_\ell^{(j)} \neq 0\big \},$$ where $\widehat{\be}^{(j)}$ is the output of lasso$(T_j(y)\sim X)$, which efficiently computes the solution of \beq \widehat{\be}^{(j)} = \operatorname*{arg\,min}_{\be^{(j)} \in \cR^{p}} \Big\{ \big\| T_j(y)\,-\,X \be^{(j)}\big\|^2_2 + \lambda \| \be^{(j)}\|_1 \Big\}. \eeq Accordingly, $\bO_1$ is the collection of location-informative variables, $\bO_2$ contains all the scale-informative ones, $\bO_3$ variables are responsible for predicting change in skewness of $Y,$ etc.

Step 3. To make the analysis doubly-robust and to allow for non-linearity of the features, simply substitute $X$ by the LP-transformed matrix $T_X$, as we did in Sec. (ref).

{\bf Example 8.} Online news popularity data, continued. Which attributes are most important for predicting the popularity of a given article? More specifically, which variables affect the changing shape of the conditional distribution? Our approach consists of four main steps:

$\bullet$ First, determine which shape parameters of $f_{Y|X=x}(y)$ are changing with $x$? In Fig. (ref), we have already seen that scale and skewness are the two principal dynamic components, in addition to location.

$\bullet$ Second, to determine the variables that impact the changing location, scale, and skewness we perform the following LP-lasso regression\footnote{see Appendix A.4.4 for a related discussion on `robust lasso.'} \[\texttt{lasso}\big(T_j(y) ~\sim~T_1(x_1) + \cdots+T_1(x_p)\big), ~~\text{for}~j=1,2,3,\] where $T_k(x_\ell):=T_k(x;\wtF_\ell)$. Store the selected features in the set $\bO_j$.

$\bullet$ Third, we summarize the findings by ranking the variables using the Ovis index; here we use $k=1$ and $j=3$ in the formula (ref). The top $10$ features are displayed in Fig. (ref), which shows the nature of contributions of different variables. For example, the variable `is_weekend' (whether the news article was published on the weekend or not) plays a dual role of being both a location and scale informant. Consequently, our multi-layered robust method provides a refined understanding of the “which and how” aspects of feature selection.

$\bullet$ Finally, one can perform targeted feature screening. An investigator can specifically query which variables are mostly responsible for heteroscedasticity, or change in symmetry, etc. Table (ref) shows the top five variables in each category: location, scale, and skewness. These insights will ultimately help the applied researchers to better understand the nature of the association between response $Y$ and $X=(X_1,\ldots,X_p)$.

The XYZ Problem: Distributional Impact Analysis

In application fields such as healthcare, economics, and social sciences, data often arise in XYZ format, where $Z$ is the binary treatment variable, $Y$ is the response variable, and $X$ is the (possibly large) collection of covariates. One such problem is discussed below.

{\bf Example 9}. Rosner's FEV data rosner95book. Table (ref) describes the data. The main interest lies in understanding the impact of smoking on the distribution of FEV as a function of age. This is important for individualized custom-tailored decision-making, where a treatment might be locally-effective for a sub-population (characterized by certain values of $x$) without being globally effective for the whole heterogeneous population.

Our model for conditional distribution of $Y$ given $X=x$ and $Z=z$ is \beq f_{Y|x,z}(y)=f_Y(y) \big\{ 1\,+\,\sum_j \LP_{j|x,z} T_j(y;F_Y)\big\}.\eeq Distributional impact analysis aims to quantify (how much) and characterize (in which ways) the differences between $f_{Y|x,1}(y)$ and $f_{Y|x,0}(y)$. To develop a measure of distributional impact, first note that from equation (ref), both $f_{Y|x,1}$ and $f_{Y|x,0}$ have the same pivot $f_Y(y)$. Thus the deviation between them directly depends on the distance between the contrast densities $d_{x,0}$ and $d_{x,1}$.

The proof is in Appendix A.3. This result motivated us to define the following measure.

We are now ready to apply this procedure on the FEV data. To estimate the LP-coefficients in (ref) we apply our UPM framework (algorithm of Section 2.3.2) with gradient-boosting machine as the base learner. We chose gbm, as it efficiently computes the interactions\footnote{ Treatment$\times$covariates interactions allow the treatment effect to vary among individuals with covariates.} between the treatment and covariates. The left panel of Fig. (ref) shows the DIF statistic for ages $12$ to $15$. The stunning fact about this graph is that it simultaneously answers which $x$-groups are most impacted by the treatment and how they are impacted. For example, at age $14$ variability is the dominant mode of difference, whereas at age $15$, location is the main source of difference. Overall, the DIF barplot does a modest job of capturing the heterogeneous impact of teenage smoking on lung function across different age groups.

Quantile Regression

Quantile regression, pioneered by koenker1978, has become a pervasive tool in a host of real-world application areas. Despite the great progress made in the last 40 years koenker2017quantile, some practical challenges still remain: how to develop a nonparametric quantile regression method that can scale to high-dimensional problems? How to efficiently incorporate nonlinearity? How to ensure that the estimated conditional quantile curves will not cross each other? Our statistical learning theory can offer some realistic solutions to these frontier problems of quantile regression. But before getting into that, let us start with a slightly broader context. Our “interlinked” statistical learning architecture (as depicted in Fig. (ref)) has two notable consequences for the practice of nonparametric data modeling: \vskip.1em $\bullet$ Our technology opens up brand-new ways of building “ML-powered” statistical models that are simultaneously flexible and scalable\footnote{By taking advantage of high-performance specialized ML-hardware and software architecture.} for large($n,p$) problems, where classical nonparametric methods are not a viable option. \vskip.1em

$\bullet$ Our UPM learning architecture is designed using a high-level universal language that is agnostic to the specific ML method, to ensure the portability and durability of the technology.\footnote{It is important to keep in mind that ML algorithms are growing at a staggeringly fast rate. In these circumstances, it is nor practical to develop `retail' strategies on a case by case basis for each ML-procedure. A practical way out is to build a unified learning ecosystem.}

\vskip.1em To perform UPM-based quantile regression, simply extract the appropriate percentiles from the conditional density estimate $\hf_{Y|X=x}(y)$. The beauty of our method is that we don't have to develop theories specialized to a certain kind of machine learning method like random forest meinshausen2006quantile,athey2019grf or neural network cannon2011quantile--it's all done using a single generic procedure; see Appendix A.4.3 and Fig. (ref) for more details.

\vskip.1em {\bf Example 10}. Dutch Boys data fredriks2000body. This dataset is a part of the Fourth Dutch Growth Study, which comprised of observations on age and BMI of $N=7294$ Dutch boys. The goal is to estimate age-specific reference growth curves based on conditional quantile function $\whQ_{Y|X=x}(u)$. Fig (ref) displays our result where we have used lasso (as implemented in the glmnet R-package) as the base learner inside UPM engine. Following the current World Health Organisation (WHO) recommendation, we have used $u=(.03, .15, .50, .85, .97)$ for constructing the cross-sectional reference growth charts.

Prediction Interval Estimator

One way to summarize the conditional density is through prediction interval (PI), which provides a concise view of the most likely values of the response variable\footnote{Knowing the range of values of $Y$ as opposed to a single-point estimate, helps decision-makers choose proper action by carefully evaluating the degree of uncertainty (length of the PIs).}. For a specified level $\al$, the goal is construct an interval that covers no less than ($1-\al$) of the probability mass of $f_{Y|X=x}(y)$. How to construct such PIs? We discuss three different procedures. The real challenge is to make it as narrow as possible, while maintaining the desired coverage.

$\bullet$ Quantile-based PI. As noted by meinshausen2006quantile, one can use quantile regression to build PIs. A $100(1-\al)$% prediction interval for $Y$ given $X=x$ can be estimated as \beq {\rm qPI}_{1-\al}(x) = \left[Q_{Y|X=x}(\al/2), Q_{Y|X=x}(1-\al/2)\right]. \eeq Naturally, the width (precision) of the quantile-PIs vary over the covariate space.

$\bullet$ Standard-error-based PI. If we are brave enough to assume some parametric form of the conditional density, then we can greatly simplify the expression for the PIs. For example, under Gaussianity assumption, one can express (ref) in the following compact form, just involving conditional mean $\mu(x)$ and standard-deviation: \beq {\rm gPI}_{1-\al}(x) = \left[ \mu(x) \pm z^{(\al/2)} \si_{|x} \right] \eeq where $\Phi^{-1}(\al)=z^{(\al)}$, $ \Var(Y|X=x)=\si^2_{|x}$, and `g' denotes the Gaussian-PI.

$\bullet$ Highest-density PI. This is the shortest possible interval with a given coverage probability boxbook73, defined as \beq {\rm hdPI}_{1-\al}(x) = \left\{y: f_{Y|x=x}(y) > \tau_{1-\al}\right\} \eeq where $\tau_{1-\al}$ is the largest number satisfying \[\int_{y \,\in\, {\rm hdPI}_{1-\al}(x)} f_{Y|X=x}(y) \dd y = 1-\al.\] Another nice property of hdPIs is that any point within the interval has a higher density than any other point outside it, which justifies their name.

Based on their operation principles, these three kinds of PIs can be broadly categorized into two groups: the vertical and the horizontal approach.

{\bf Example 11} Butterfly data, continued. All three prediction intervals are shown in Fig. (ref)(a). Three main conclusions: (i) The shortest among all is the highest-density PI, which take the form of disjoint subintervals--one for each local mode. (ii) The Gaussian assumption of {\rm gPI} seems too restrictive for real-world data analysis. (iii) Both {\rm gPI} and {\rm qPI} include a low-density valley area around $x=0$ at the cost of other, more-likely values around $x=\pm 2$.

{\bf Example 12}. The Auto-MPG data. The task is to predict the miles-per-gallon (MPG) gas consumption of an automobile based on $p=7$ features including horsepower, weight, and acceleration. The dataset was used in the 1983 American Statistical Association Exposition and is available in the UCI Machine Learning Repository. Here, we are mainly interested in predicting the MPG of a specific instance--row number $390$, which corresponds to the 1982 Dodge Rampage car. The UPM-predicted uncertainty distribution is displayed in Fig. (ref)(b). The bimodality of the estimated conditional density reflects the fact that 1982 Dodge Rampage is a hybrid two-door car--a mix between a passenger car and a truck. The next notable thing is the length of the PIs: The standard error-based gPI is almost three times wider than the hdPI, and thus provides an overly pessimistic assessment of an accurate prediction. This also sheds some light on the apparently paradoxical fact (as noticed by wager2014confidence) why bootstrap standard error-based Gaussian confidence interval yields a large error bar for the Dodge Rampage example, even when the prediction was dead on.

Modern Applied Statistics: Theory, Practice, and Pedagogy

A theoretical research, which has no relevance to practice and teaching, is not persuasive enough to be taken as fundamental. To us, it is very important to know whether our research was able to link (at least partially) the triad of theory-practice-teaching. In this article, through several examples, we described `the joy of systematic data analysis' derived from a general theory, which has the following implications for the practice and pedagogy of statistics:

$\bullet$ Our theory makes it easy: to apply (arise from simplification) and see the connections (arise from unification) between different statistical methods.

$\bullet$ Our theory provides a holistic training\footnote{For more discussion on this topic see Deep17Teaching.} by introducing a large variety of statistical topics (embracing different cultures: algorithmic, parametric, nonparametric, information-theoretic, robust, exploratory) in an organized manner through a common semantics. This is extremely important to produce broad integrator foxes rather than highly specialized hedgehogs.\footnote{Readers may find it worthwhile to read Philip Tetlock's book on `Why Foxes Are Better Forecasters Than Hedgehogs,' which is based on a 20-year study with 284 experts from diverse fields, including government officials, professors, journalists, and others. As Richard Feynman said “Science is not a specialist business.”} Ultimately, what matters is whether our theory provides a faster and easier (lazier?) way to learn the fundamentals of statistics. We believe--yes, it does.\footnote{At least it has a better chance than any other currently known global theory of data analysis. We derive confidence from the ability of our theory to tackle statistical problems as diverse as time-series analysis D12e, copula modeling D20copula, empirical Bayes deep18nature, graph-theory Deep19SGT, multiple testing Deep18DMT, deep16LSSD, high-dimensional data analysis DeepLPKsample19, large-scale distributed learning deep16MetaLP, etc. This is an ongoing and growing movement with the mission to structure the field in an understandable and effective way.}

70 Years of Statistical Machine Learning

Breiman's 2001 paper was highly influential and hailed as a heroic effort to challenge the status quo in a “David and Goliath” style, where two diametrically opposite modeling cultures--machine learning and traditional parametric statistics--battled for the throne. However, it is much less well-known among practitioners (including the current generation of deep-learning engineers) that machine learning has strong historical roots in (nonparametric) Statistics. In an unpublished US Air Force School of Aviation Medicine report in 1951, two UC Berkeley Statisticians (students of Jerzy Neyman) Evelyn Fix and J.L. Hodges, Jr., introduced a non-parametric method for pattern classification, which is now called the k-nearest neighbor (knn) method. This was a groundbreaking paper that marked the beginning of statistical machine learning\footnote{Interestingly, just after that, Arthur Samuel came up with the phrase “Machine Learning” in 1952.}. Fix and Hodges's work was a precursor to the celebrated kernel-based methods and decision tree-based algorithms, which revolutionized the world of machine learning. The time has come that we bring machine learning back to its statistical roots.\footnote{The next revolution of machine learning will need both power of computing and power of statistical thinking (and little bit of common sense).} \nocite{fix51}

\vskip.25em

The current paper challenges Breiman's “two-culture” model in favor of an integrated learning system where different cultures can be “glued” together by some underlying deeper principles.\footnote{There is no doubt that in the coming years we will witness many more new varieties and cultures of data analysis, but in the midst of those temporary excitements, we should not forget to keep the discipline internally consistent. Diversity without unity would be a complete mess --not healthy for our profession.} Through numerous examples, we have shown that the two cultures are not contradictory but complementary to each other, which when combined properly, can yield a more powerful learning technology that is simultaneously flexible, scalable, and explainable.

\vskip.25em

If there is one thing that our work has made apparent is: there's plenty of room in the middle--tremendous opportunities lie at the boundaries between the two cultures. In this paper we have described a new theory for constructing a universal interface between statistics and machine learning. What we have proposed is by no means the only possible or even the “best” theory; it is merely a provisional one. As we go forward, it will be very important to come up with new innovative theories that can reconcile and build new links between various cultures of statistical modeling. All of this is still in its nascent stages, but definitely poised for rapid progress in the coming decade.

Acknowledgement

This paper celebrates the 200th anniversary of Gauss-Markov least-square regression and 70th anniversary of statistical machine learning\footnote{These two historical moments symbolize the birth of `two cultures' of data modeling--linear parametric and flexible algorithmic modeling cultures.}, by putting forward a modern unified perspective.

{1.6ex}

\addcontentsline{toc}{section}{6 Appendix}

A.1. Software

We have developed an R Software LPMachineLearning\footnote{A preliminary version is available from the authors. A more upgraded one is currently under development.} to perform all the tasks outlined in the paper. We hope this software will encourage applied data scientists to apply our method for their real prediction problems.

\vskip.6em

A.2. LP-Orthonormal System: Construction and Properties

LP-basis functions are specially-designed orthonormal system of rank-polynomials. They are inherently `nonparametric,' constructed in a fully data-driven way, not pre-defined like classical polynomials. For a mixed (either discrete or continuous) random variable $Z$ with distribution function $F_Z$, we describe the construction of LP-basis functions. \vskip.4em

{\bf Construction Algorithm}. A fully-automated construction of LP-orthonormal rank-polynomials of $Z$ consists of the following steps: \vskip.3em Step 1. Mid-Distribution Transform: The mid-distribution function of $Z$ is defined as \beq \Fm_Z(z)=F_Z(Z)-\frac{1}{2}p_Z(z), \eeq where $p_Z(z)$ is the probability mass function. \vskip.2em Step 2. Standardizing $\Fm_Z(Z)$: The random variable $\Fm_Z(Z)$ has mean $\Ex[\Fm_Z(Z)]=.5$ and variance $\Var[\Fm_Z(Z)]=\frac{1}{12}\big( 1- \sum_z p_Z^3(z) \big)$. Define the 1st-order LP-basis function as \beq T_1(z;F_Z) = \dfrac{\sqrt{12}\big\{\Fm_Z(z) - 1/2\big\}}{\sqrt{1-\sum_z p_Z^3(z)}}, \eeq \vskip.1em Step 3. Gram-Schmidt Orthonormalization: Construct high-order polynomials by applying Gram-Schmidt orthonormalization (see Appendix A of D20copula for a quick refresher on Gram-Schmidt process) on the set of functions of the power of $T_1(Z;F_Z)$.

\vskip.4em {\bf Some properties}. Here we list some important properties of the They satisfy the following properties of the LP-orthonormal systems:

1) LP-bases are orthonormal with respect to the measure (weighting function) $F_Z$: \[ \int T_j(z;F_Z) \dd F_Z =0,\,~~\text{and} \,\,\int T_j(z;F_Z) T_k(z;F_Z) \dd F_Z =\delta_{jk}.\] For related discussions, see Remark 4 of D20copula1.

\vskip.35em 2) For real-data analysis, construct empirical LP basis (in short eLP basis) $\{T_j(z;\wtF_Z)\}_{j=1,2\ldots,m}$, where $m$ is strictly less than the number of unique values in the sample $\{z_1,\ldots, z_n\}$. eLP bases are orthonormal system of functions with respect to the discrete empirical measure $\wtF_Z$.

\vskip.15em

3) For $Z$ continuous, it is not difficult to show that the LP-bases reduces to the following universal shape: \beq T_j(z;F_Z) := \Leg_j \circ \,F_Z(z),\eeq where $\Leg_j(u)$ denotes $j$th shifted orthonormal Legendre polynomials on $[0,1]$ and `$\circ$' denotes composition of functions.\footnote{In our regression context the response variable $Y$ is continuous, thus admits this analytic LP-basis form. However, this one-to-one correspondence between LP-polynomials and Legendre polynomials of rank-transform (also know as probability integral transform) $F_Z(Z)$ is not true for discrete $Z$, i.e., for the covariates $X$.} Few top LP-polynomials for $Z$ continuous are given below: \bea T_1(z;F_Z)&=&\Leg_1\big( F_Z(Z) \big)\,=\, \sqrt{12}\big( F_Z(z) - 1/2 \big) \\ T_2(z;F_Z)&=&\Leg_2\big( F_Z(Z) \big)\,=\,\sqrt{5}\big( 6 F^2_Z(Z) - 6 F_Z(Z) + 1 \big) \\ T_3(z;F_Z)&=&\Leg_3\big( F_Z(Z) \big)\,=\,\sqrt{7}\big( 20 F^3_Z(Z) - 30 F^2_Z(Z) + 12 F_Z(Z) -1 \big). \eea and so on. For more information, ses DD20181,Deep17LPMode1. This will be useful for constructing LP-polynomials for any general pivot density $f_0$ of the response variable. \vskip.15em

4) Derivation of empirical $T_1(z;\wtF_Z)$ from a sample of all unique $z_1,\ldots,z_N$ (i.e, realization from a continuous $Z$). In this case, the empirical probability mass function (pmf) is $\tp_Z(z_i)=1/N$, and the empirical cdf is $\wtF_Z(z_i)=R_i/N$, where $R_i$=rank($z_i$). This immediately implies (following eq. (ref)) \beq T_1(z;\wtF_Z) = \sqrt{\frac{12}{N^2-1}}\Big(R_i - \frac{N+1}{2}\Big). \eeq Since, \[ \tFm_Z(z_i)= \wtF_Z(z_i)- \frac{1}{2}\tp_Z(z_i)= \frac{1}{N} (R_i - 1/2).\] and the correction factor \[ 1-\sum_i \tp_Z^3(z_i) = 1 - \frac{1}{N^2}. \] Eq. (ref) makes it clear why we call LP-bases are polynomials of ranks.

\vskip.15em 5) The shapes of LP-polynomials are data-adaptive. The left panel of Fig. (ref) displays the top four eLP-basis functions for the covariate weight, taken from the Auto-MPG data (Example 12 in the main paper). The weight variable has very little ties, which is the reason why the shapes of LP-bases completely match with eq. ((ref)-(ref)). Contrast this with the right panel of Fig. (ref), which shows top four eLP bases for the variable acceleration. Due to the presence of a large number of ties, it takes a unique piecewise-constant shape.

\vskip1em

A.3. Remaining Proofs

\vskip.35em {\bf Theorem (ref)}: Since the conditional entropy $H(Y|X)=\int H(Y|X=x) \dd F_{X}(x)$, we have \beas H(Y) - H(Y|X) &=& - \int f_Y(y) \log f_Y(y) \dd y + \iint f_{Y|x}(y) \log f_{Y|x}(y) \dd y \dd F_X(x)\\ &=& \int \Big \{ \int f_{Y|x}(y) \log \frac{f_{Y|x}(y)}{f_Y(y)} \dd y \Big\} \dd F_X(x)\\ &=& \Ex_X{\rm KL}(d_X;\mathbbm{1}_{[0,1]}). \eeas The last equality follows from Theorem (ref). \vskip.35em {\bf Theorem (ref)}: The Kruskal–Wallis statistic (ref) is given by (up to a negligible factor $(n-1)/n$) \beq KW = \sum_{\ell=1}^k n_\ell \Bigg\{ \sqrt{\frac{12}{N^2-1}}\big( n_\ell^{-1}\sum_{i\in G_\ell} R_i - \frac{N+1}{2} \big) \Bigg\}^2 \eeq Now note that \beq \tLP_{1|\ell} = \Ex[T_1(Y;\wtF_Y)|X=\ell] = \frac{1}{n_\ell} \sum_{i \in G_\ell} \sqrt{\frac{12}{N^2-1}}\Big(R_i - \frac{N+1}{2}\Big). \eeq The first equality follows from (ref) and the second one from (ref). Combining (ref) and (ref) we get the first desired result \beq KW = \sum_{\ell=1}^k n_\ell \,\big| \tLP_{1|\ell} \big|^2.\eeq Next, to show that (ref) is actually equals to $NR_1^2$, note that \[R_1^2 ~=~\frac{\Var(E[T_1(Y;F_Y)|X])}{\Var(T_1(Y;F_Y))},\] where the denominator is $1$ by the construction of LP-basis functions, and the numerator is $\sum_{\ell=1}^k \pi_\ell \,\big| \tLP_{1|\ell} \big|^2$, $\pi_\ell$ is $\Pr(X=\ell)=n_\ell/N$. This is because \[\Ex[\LP_{1|X}]~=~\Ex\left[\Ex[T_1(Y;F_Y)|X]\right]~=~\Ex[T_1(Y;F_Y)]=0.\]

{\bf Theorem (ref)}: Define the LP unit-basis functions as \beq S_j(u;F_Y)=T_j(Q_Y(u);F_Y) {\rm for} 0<u<1,\eeq which satisfy the following orthonormality relations (immediate from eq. (ref)) \beq \int_0^1 S_j(u;F_Y) S_k(u;F_Y) \dd u\,=\,\delta_{jk}.\eeq

Substitute $F_Y(y)=u$ in (ref) to express the contrast density function in the quantile domain: \[d_{x,z}(u)~=~1+\sum_{j} \LP_{j|x,z} S_j(u;F_Y),~0<u<1.\] It is now easy to see that the $L^2$ distance between $d_{x,0}$ and $d_{x,1}$ can be written as \[\big \| d_{x,0} - d_{x,1} \big\|_2^2 \,=\, \sum_{j} \big| \LP_{j|x,0}\,-\, \LP_{j|x,1}\big|^2 \] due to (ref). Hence proved. \vskip1em

A.4. Additional Remarks

\vskip.35em {\bf Remark A.4.1} Koenker's Quantile regression for Butterfly data: Fig. (ref) shows the quantile regression curves for the butterfly data, which implements the algorithm proposed in koenker1987algorithm.

\vskip.35em {\bf Remark A.4.2} Smoothness of Regression via Rank-transform. Fig. (ref) shows house prices (median value of in 1000s) in Boston metropolitan area as function of crime (per capita crime rate). The data were collected in 1970 and available in the UCI Machine Learning Repository. The left plot shows the data in the xy-domain, and the right plot shows the same data in the quantile domain where we replaced $x$ by its rank-transform $\wtF_X(x)$. Red curves are the spline smoother (estimated using smooth.spline() R-function). Few interesting points to note:

$\bullet$ Smoothness: Regression via rank-transform (see Theorem (ref)) seems to produce much smoother (parsimonious) regression curve compared to the rough zigzag xy-domain estimate. Is it universal that nature reveals herself in a more parsimonious way in the quantile domain? Apart from intuitive understanding, we still don't know any concrete mathematical explanation for this surprising phenomenon.

$\bullet$ Data-sparsity: The quantile-domain treatment seems to address the `data sparsity' (that exists for $x>30$) problem quite well. Surprisingly, after quantile-transformation, a data-sparse regression problem turns into a data-dense one.

$\bullet$ Visualization: The rank-transform scatter makes it far easy to visualize and understand the relationship between $x$ and $y$.

$\bullet$ Robust+Smooth: This shows an additional benefit of X-robustness (Theorem (ref)) as a tool for promoting smoothness.

\vskip.35em {\bf Remark A.4.3} ML-assisted Quantile regression: comparison. The first row of Fig. (ref) shows the generalized random forests-based athey2019grf1 and gradient boosted quantile regression estimates. They are specially-designed ML-algorithms to produce the desired conditional quantile curves. The second row shows our UPM-based versions, which are derived from a general scheme described in Section (ref) of the main paper.

\vskip.35em {\bf Remark A.4.4} Robust Lasso: As already noted, the online news data contain an outlier--the news article in row \#31,038. How does it impact traditional estimation and feature selection methods? As we will see, a single data point can have a devastating impact on the overall modeling. Fig. (ref) (left panel) shows the lasso estimated $\widehat{\be}$ under two scenarios: first, on the outlier-removed data, and second, on the full data (with outlier). Notice the estimated coefficient values (as well as their sparsity-levels) are dramatically affected because of just one data point. One way to address this fragility is to apply lasso on the first-order LP-transformed $T_X=[T_1(x;\wtF_{X_1}),\ldots, T_1(x;\wtF_{X_p})]$, instead of raw feature matrix $X$. The result is displayed on the right panel, which shows extraordinary consistency with (almost) no impact on the estimated coefficients. This mid-rank transform based robust lasso method could be a reliable alternative to deal with noisy messy data.