Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
97,491 characters · 12 sections · 47 citation commands
Ordered Correlation Forest
\doublespacing
Categorical outcomes with a natural order, often referred to as ordered non-numeric outcomes, are commonly observed in empirical studies across the social sciences. For example, happiness research typically employs large surveys to collect self-evaluations of subjective well-being frey2002can, and health economics is heavily based on self-assessments in several health domains peracchi2012heterogeneity, peracchi2013heterogeneous. These outcomes are usually measured on a discrete scale with five or ten classes, where the classes can be arranged in a natural order without any knowledge about their relative magnitude.
Ordered choice models, such as ordered logit and ordered probit, are frequently used to analyze the relationship between an ordered outcome and a set of covariates greene2010modeling. These models target the estimation of the conditional choice probabilities, which represent the probability that the outcome belongs to a certain class given the values of the covariates. However, they are limited by their dependence on parametric and distributional assumptions that are often based on analytical convenience rather than knowledge about the underlying data generating process. As a result, econometricians may need to consider alternative techniques to produce more accurate and reliable predictions.
This paper introduces a novel machine learning estimator specifically optimized for handling ordered non-numeric outcomes. Employing traditional machine learning estimators “ off-the-shelf" can result in biased and inefficient estimation of conditional probabilities. This is because classification algorithms do not leverage the ordering information embedded in the structure of the outcome, and regression algorithms treat the outcome as if it is measured on a metric scale.\footnote{\ For comprehensive overviews of traditional classification and regression algorithms, the reader is referred to hastie2009elements and efron_hastie_2016.} The proposed estimator is designed to mitigate the biases that traditional methods can introduce, ultimately resulting in enhanced predictive performance.
The proposed estimator, named the ordered correlation forest, adapts a standard random forest splitting criterion breiman2001random to the mean squared error relevant to the specific estimation problem at hand. The new splitting rule is then used to build a collection of forests, each estimating the conditional probability of a single class. After constructing the individual trees within each forest, the ordered correlation forest employs an unbiased estimator of conditional probabilities within each leaf. Model consistency is ensured, as the predictions always fall within the unit interval by construction. To estimate the covariates' marginal effects, the ordered correlation forest utilizes a nonparametric approximation of derivatives lechner2019random.
Under an “ honesty" condition athey2016recursive, the ordered correlation forest inherits the asymptotic properties of random forests, namely the consistency and asymptotic normality of their predictions wager2018estimation. Honesty is a subsample-splitting technique that requires that different observations are used to place the splits and compute leaf predictions and is crucial to achieving consistency of the random forest predictions.
The particular honesty implementation used by the ordered correlation forest allows for a weight-based estimation of the variance of the predicted probabilities. This is achieved by rewriting the random forest predictions as a weighted average of the outcomes athey2019generalized. The weights, which are obtained for the predicted probabilities, can be properly transformed to obtain standard errors for the covariates' marginal effects lechner2022modified, lechner2019random. We can then use the estimated standard errors to conduct valid inference about the marginal effects as usual, e.g., by constructing conventional confidence intervals.
The rest of the paper unfolds as follows. Section (ref) provides a brief overview of ordered choice models and discusses some alternative estimation strategies. Section (ref) presents the ordered correlation forest, explaining estimation and inference about the statistical targets of interest. Section (ref) uses synthetic data to compare the ordered correlation forest with alternative estimators and evaluate its performance in estimating and making inference about the covariates' marginal effects. Section (ref) provides further comparisons with alternative estimators using real data. Section (ref) concludes.
Ordered choice models are a class of statistical models used to analyze the relationship between an ordered non-numeric outcome $Y_i$ and a set of covariates $W_i$ mccullagh1980regression. These models are typically motivated by postulating the existence of a latent and continuous outcome variable of interest $Y_i^*$, assumed to obey the following regression model peracchi2014econometric:
where $W_i$ consists of a set of raw covariates, $g \left( \cdot \right)$ is a potentially non-linear regression function, and $\epsilon_i$ is independent of $W_i$ and has cumulative distribution $F \left( \cdot \right)$. Then, an observational rule links the observed outcome $Y_i$ to the latent outcome $Y_{i}^*$ using unknown threshold parameters $- \infty = \zeta_0 < \zeta_1 < \dots < \zeta_{M - 1} < \zeta_M = \infty$ that define intervals on the support of $Y_i^*$, with each interval corresponding to one of the $M$ categories or classes of $Y_i$:
Although the $M$ classes have a natural ordering, they are not measured on a cardinal scale. This limits our ability to make precise quantitative comparisons.
Researchers are typically interested in the estimation of the conditional choice probabilities, defined as:
However, the marginal effect of the $j$-th covariate on $p_m \left( \cdot \right)$ is a more interpretable measure for ordered choice models. The marginal effect is defined differently depending on whether the $j$-th covariate is continuous or discrete:
where $w_j$ is the $j$-th element of the vector $w$ and $\lceil w_j \rceil$ and $\lfloor w_j \rfloor$ correspond to $w$ with its $j$-th element rounded up and down to the closest integer. We can summarize the marginal effects in various ways, such as computing the marginal effect at the mean $\nabla^j p_m \left( \bar{w} \right)$, with $\bar{w}$ denoting a vector of means. Alternatively, we can compute the marginal effect at the median, the mean marginal effect, and the median marginal effect.
From ((ref)) and ((ref)), the conditional choice probabilities write as:
If the regression function $g \left( \cdot \right)$ and the distribution $F \left( \cdot \right)$ of the error term $\epsilon_i$ are known, we can estimate ((ref)) directly using standard maximum likelihood methods.
However, in many practical applications, precise knowledge of $g \left( \cdot \right)$ is not available. Instead, a common approach is to approximate it using a linear-in-parameter model belloni2011high:
where $X_i = h \left( W_i \right)$ is a $k$-dimensional vector of constructed covariates (generally the raw covariates $W_i$ plus interactions and polynomials thereof) and $V_{i, k}$ is an approximation error that is assumed to be independent of $X_i$. Substituting the linear approximation ((ref)) into ((ref)) gives:
where the random error $U_i = \epsilon_i + V_{i, k}$ depends on $k$ through the approximation error $V_{i, k}$ and has cumulative distribution $G \left( \cdot \right)$. Then, we can approximate the statistical target $p_m \left( \cdot \right)$ as follows:\footnote{\ The ultimate target of estimation is $p_m \left( \cdot \right)$. $p_m^* \left( \cdot \right)$ serves as an approximation that allows us to tackle the estimation problem as if it were parametric.}
We can impose assumptions on the distribution $G \left( \cdot \right)$ of the random error $U_i$ to estimate ((ref)) using standard maximum likelihood methods. Popular choices are the standard normal and the standard logistic distribution functions, producing the ordered probit and ordered logit models, respectively. In scenarios where $k > n$, regularization techniques such as L1- or L2-type penalization are needed.
Although easy to interpret and computationally simple, this approach features several limitations. First, it imposes strong distributional assumptions generally derived from analytical convenience rather than knowledge about the underlying data generating process. Second, it requires the specification of a linear-in-parameter model such as ((ref)) to account for non-linearities in $g(\cdot)$. Third, the estimated marginal effects have the restrictive property of single-crossing, meaning that they can change sign only once when moving from the smallest class to the largest.
Recent developments in statistical learning hastie2009elements, efron_hastie_2016 offer ways to overcome these limitations. For instance, random forest algorithms breiman2001random offer a nonparametric estimation approach that does not assume a specific error term distribution and can naturally handle non-linearities in $g \left( \cdot \right)$ without requiring a linear-in-parameter model. However, classification forests do not leverage the ordering information embedded in the structure of the outcome, and regression forests treat the outcome as if it is measured on a metric scale. Consequently, applying these algorithms “ off-the-shelf" can result in biased and inefficient estimation of conditional probabilities.
To overcome these limitations, one approach is to transform ordered non-numeric outcomes into a metric scale using scores based on the classes of the observed outcome, thus allowing us to use any regression algorithm on the transformed outcome. For example, hothorn2006unbiased propose using the midpoint values of the intervals defined on the support of the latent outcome as score values. In the cases where $Y_i^*$ is not observed, this translates into setting the scores equal to the class labels of $Y_i$. However, this assumes that the intervals are of equal length, which may not be accurate in practice. To address this issue, hornung2020ordinal proposes the ordinal forest estimator, which optimizes the class intervals and uses score values corresponding to these optimized intervals in a standard regression forest. The optimization process involves growing multiple forests using randomly generated candidate score sets, and constructing the final score values by summarizing the score sets with the smallest out-of-bag error. hornung2020ordinal shows that the ordinal forest estimator outperforms a standard regression forest that uses class labels as score values using both real and synthetic data. However, the optimization process can be computationally expensive, which may limit its practical use for large data sets or real-time applications.
Another approach involves expressing conditional probabilities as conditional expectations of binary variables, which can be estimated by any regression algorithm. One first strategy, which we label multinomial machine learning, is to express conditional probabilities as follows:
This allows us to estimate each $p_m \left( \cdot \right)$ separately by regressing the binary variable $\vmathbb{1} \left( Y_i = m \right)$ on $W_i$ using any nonparametric estimator:
Alternatively, we can specify a linear-in-parameter model to estimate the approximate target $p_m^* \left( \cdot \right)$ through parametric regression of the binary variable $\vmathbb{1} \left( Y_i = m \right)$ on $X_i$.
However, $\hat{p}_m^{MML} \left( \cdot \right)$ does not leverage the information embedded in the ordered structure of the outcome. To overcome this limitation, an alternative strategy that we label ordered machine learning expresses conditional choice probabilities as the difference between the cumulative probabilities of two adjacent classes:
with $\mu_m \left( W_i \right) := \operatorname{\mathbb{E}} \left[ \vmathbb{1} \left( Y_i \leq m \right) | W_i \right]$. Then we can estimate each $\mu_m \left( \cdot \right)$ separately by regressing the binary variable $\vmathbb{1} \left( Y_i \leq m \right)$ on $W_i$ using any nonparametric estimator and pick the difference between the cumulative probabilities of two adjacent classes to estimate $p_m \left( \cdot \right)$:\footnote{\ lechner2019random combine ordered machine learning with random forests breiman2001random and discuss how to estimate and conduct inference about marginal effects.}
As before, we can alternatively specify a linear-in-parameter model to estimate the approximate target $p_m^* \left( \cdot \right)$ through parametric regressions of the binary variables $\vmathbb{1} \left( Y_i \leq m \right)$ and $\vmathbb{1} \left( Y_i \leq m - 1 \right)$ on $X_i$.
However, $\hat{p}_m^{OML} \left( \cdot \right)$ can potentially produce negative predictions, thereby contradicting the definition of probabilities. Although we might resolve this issue by setting negative predictions to zero, such a solution is suboptimal, and an alternative estimator that does not require truncation may perform better. This paper introduces a novel estimator that leverages the ordered structure of the outcome and produces predictions that always fall within the unit interval, thus resulting in enhanced predictive performance compared to existing methods.
In this section, I discuss the implementation of the ordered correlation forest (OCF) estimator. First, I illustrate the estimation of conditional choice probabilities and marginal effects. Second, I discuss the conditions required for the consistency and asymptotic normality of OCF predictions. Finally, I show how to conduct approximate inference about the statistical targets of interest.
Similar to the ordered machine learning approach, OCF computes the prediction of conditional choice probabilities as the difference between the cumulative probabilities of two adjacent classes (see equation (ref)). However, instead of estimating $\mu_m \left( \cdot \right)$ and $\mu_{m - 1} \left( \cdot \right)$ separately, OCF internally performs this computation in a single random forest. This allows us to tie the estimation of $\mu_m \left( \cdot \right)$ and $\mu_{m - 1} \left( \cdot \right)$ to correlate the errors made in estimating these two expectations. Additionally, it avoids negative predictions.
To see the importance of correlating the estimation errors, we can decompose the mean squared error of a prediction $\hat{p}_m^{OML} \left( \cdot \right)$ at $w$ as follows:\footnote{\ This decomposition can be applied to any estimation strategy that involves calculating the difference between two surfaces. For example, lechner2022modified leverage this decomposition to estimate heterogeneous causal effects under a selection-on-observables assumption}.
where the last term is the error correlation and captures the degree to which the errors made in estimating $\mu_m \left( \cdot \right)$ and $\mu_{m - 1} \left( \cdot \right)$ are correlated:
Equation ((ref)) shows that $\hat{p}_m^{OML} \left( \cdot \right)$ is a suboptimal estimator. Besides potentially leading to negative predictions, estimating $\mu_m \left( \cdot \right)$ and $\mu_{m - 1} \left( \cdot \right)$ separately minimizes only the mean squared error terms and ignores the error correlation. Tying the estimation of $\mu_m \left( \cdot \right)$ and $\mu_{m - 1} \left( \cdot \right)$ to correlate the errors could improve estimation performance since errors that move in the same direction cancel out when taking the difference $\hat{\mu}_m \left( \cdot \right) - \hat{\mu}_{m - 1} \left( \cdot \right)$.
To address this limitation, OCF constructs a collection of forests, one for each of the $M$ classes of $Y_i$. However, rather than the standard criterion breiman2001random, OCF uses equation ((ref)) as the splitting rule to build the individual trees in the $m$-th forest. This allows the estimator to account for the error correlation that $\hat{p}_m^{OML} \left( \cdot \right)$ ignores. Intuitively, during the tree-building process, OCF anticipates that the predictions in the final leaves will involve the difference between two estimated functions. Consequently, it seeks splits that not only yield accurate estimates of $\mu_m \left( \cdot \right)$ and $\mu_{m - 1} \left( \cdot \right)$ but also take into account the correlation between the errors made in estimating these expectations.
To use ((ref)) as the splitting rule, we need to estimate its components. This, in turn, requires an estimator of $\mu_m \left( \cdot \right)$ in each node. An unbiased estimator of $\mu_m \left( \cdot \right)$ in a child node $\mathcal{C}_j \subset \mathcal{W}$ consists of the proportion of observations in $\mathcal{C}_j$ whose outcome is not greater than $m$:
This leads us to estimating $\text{MSE} \left( \check{\mu}_m \left( \cdot \right) \right)$ and $\text{EC} \left( \check{\mu}_m \left( \cdot \right), \check{\mu}_{m - 1} \left( \cdot \right) \right)$ in each node by their sample analogs:
Then, in the $m$-th forest, OCF constructs individual trees by recursively partitioning each parent node $\mathcal{P} \subseteq \mathcal{W}$ into two child nodes $\mathcal{C}_1, \mathcal{C}_2 \subset \mathcal{P}$ such that the following minimization problem is solved:
Once the recursive partitioning stops, each tree in the $m$-th forest unbiasedly estimates $p_m \left( \cdot \right)$ at $w$ by computing the proportion of observations in the same leaf as $w$ whose outcome equals $m$:
where $L_{m, b} \left( w \right)$ is the set of observations falling in the same leaf of the $b$-th tree as the prediction point $w$. The predictions from each tree are then averaged to obtain the forest predictions:\footnote{\ It may be necessary to perform a normalization step to ensure that $\sum_{m = 1}^M \hat{p}_m^{OCF} \left( w \right) = 1$. This is true also for $\hat{p}_m^{MML} \left( \cdot \right)$ and $\hat{p}_m^{OML} \left( \cdot \right)$.}
where $b = 1, \dots, B_m$ indexes the trees in the $m$-th forest. In contrast to ordered machine learning, OCF ensures model consistency, as the predictions $\hat{p}_m^{OCF} \left( \cdot \right)$ always fall within the unit interval by construction.
Estimation of marginal effects proceeds as proposed by lechner2019random. For discrete covariates, we can plug an estimate $\hat{p}_m^{OCF} \left( \cdot \right)$ of $p_m \left( \cdot \right)$ into equation ((ref)) to have a straightforward estimator of $\nabla^j p_m \left( \cdot \right)$:
For continuous covariates, we use a nonparametric approximation of the infinitesimal change in $w_j$:
where $\widehat{\lceil w_j \rceil}$ and $\widehat{\lfloor w_j \rfloor}$ correspond to $w$ with its $j$-th element set to $\widebar{w}_j = w_j + \omega \, \sigma_j$ and $\@ifnextchar^{{\if@single{w}{ \begingroup \def\mathaccent#w#0{ \let\mathaccent\save@mathaccent \if12 \let\macc@nucleus \fi \setbox\z@\hbox{$\macc@style{\macc@nucleus}_$} \setbox\tw@\hbox{$\macc@style{\macc@nucleus}_$} \dimen@\wd\tw@ \advance\dimen@-\wd\z@ \divide\dimen@ 3 \@tempdima\wd\tw@ \advance\@tempdima-\scriptspace \divide\@tempdima 10 \advance\dimen@-\@tempdima \ifdim\dimen@>\z@ \dimen@0pt\fi \rel@kern{0.6}\kern-\dimen@ \if11 \underbar{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@} \advance\[email removed]\dimexpr\macc@kerna \let\final@kern0 \ifdim\dimen@<\z@ \let\final@kern1\fi \if\final@kern1 \kern-\dimen@\fi \else \underbar{\rel@kern{-0.6}\kern\dimen@w} \fi } \macc@depth\@ne \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax} \macc@set@skewchar\relax \let\mathaccentV\macc@nested@a \if11 \macc@nested@a\relax111{w} \else \def\gobble@till@marker#w\endmarker{} \futurelet\gobble@till@markerw\endmarker \ifcat\noexpand A\else \def{} \fi \macc@nested@a\relax111{} \fi \endgroup }{ \begingroup \def\mathaccent#w#0{ \let\mathaccent\save@mathaccent \if22 \let\macc@nucleus \fi \setbox\z@\hbox{$\macc@style{\macc@nucleus}_$} \setbox\tw@\hbox{$\macc@style{\macc@nucleus}_$} \dimen@\wd\tw@ \advance\dimen@-\wd\z@ \divide\dimen@ 3 \@tempdima\wd\tw@ \advance\@tempdima-\scriptspace \divide\@tempdima 10 \advance\dimen@-\@tempdima \ifdim\dimen@>\z@ \dimen@0pt\fi \rel@kern{0.6}\kern-\dimen@ \if21 \underbar{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@} \advance\[email removed]\dimexpr\macc@kerna \let\final@kern0 \ifdim\dimen@<\z@ \let\final@kern1\fi \if\final@kern1 \kern-\dimen@\fi \else \underbar{\rel@kern{-0.6}\kern\dimen@w} \fi } \macc@depth\@ne \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax} \macc@set@skewchar\relax \let\mathaccentV\macc@nested@a \if21 \macc@nested@a\relax111{w} \else \def\gobble@till@marker#w\endmarker{} \futurelet\gobble@till@markerw\endmarker \ifcat\noexpand A\else \def{} \fi \macc@nested@a\relax111{} \fi \endgroup }}}{\if@single{w}{ \begingroup \def\mathaccent#w#1{ \let\mathaccent\save@mathaccent \if12 \let\macc@nucleus \fi \setbox\z@\hbox{$\macc@style{\macc@nucleus}_$} \setbox\tw@\hbox{$\macc@style{\macc@nucleus}_$} \dimen@\wd\tw@ \advance\dimen@-\wd\z@ \divide\dimen@ 3 \@tempdima\wd\tw@ \advance\@tempdima-\scriptspace \divide\@tempdima 10 \advance\dimen@-\@tempdima \ifdim\dimen@>\z@ \dimen@0pt\fi \rel@kern{0.6}\kern-\dimen@ \if11 \underbar{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@} \advance\[email removed]\dimexpr\macc@kerna \let\final@kern1 \ifdim\dimen@<\z@ \let\final@kern1\fi \if\final@kern1 \kern-\dimen@\fi \else \underbar{\rel@kern{-0.6}\kern\dimen@w} \fi } \macc@depth\@ne \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax} \macc@set@skewchar\relax \let\mathaccentV\macc@nested@a \if11 \macc@nested@a\relax111{w} \else \def\gobble@till@marker#w\endmarker{} \futurelet\gobble@till@markerw\endmarker \ifcat\noexpand A\else \def{} \fi \macc@nested@a\relax111{} \fi \endgroup }{ \begingroup \def\mathaccent#w#1{ \let\mathaccent\save@mathaccent \if22 \let\macc@nucleus \fi \setbox\z@\hbox{$\macc@style{\macc@nucleus}_$} \setbox\tw@\hbox{$\macc@style{\macc@nucleus}_$} \dimen@\wd\tw@ \advance\dimen@-\wd\z@ \divide\dimen@ 3 \@tempdima\wd\tw@ \advance\@tempdima-\scriptspace \divide\@tempdima 10 \advance\dimen@-\@tempdima \ifdim\dimen@>\z@ \dimen@0pt\fi \rel@kern{0.6}\kern-\dimen@ \if21 \underbar{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@} \advance\[email removed]\dimexpr\macc@kerna \let\final@kern1 \ifdim\dimen@<\z@ \let\final@kern1\fi \if\final@kern1 \kern-\dimen@\fi \else \underbar{\rel@kern{-0.6}\kern\dimen@w} \fi } \macc@depth\@ne \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax} \macc@set@skewchar\relax \let\mathaccentV\macc@nested@a \if21 \macc@nested@a\relax111{w} \else \def\gobble@till@marker#w\endmarker{} \futurelet\gobble@till@markerw\endmarker \ifcat\noexpand A\else \def{} \fi \macc@nested@a\relax111{} \fi \endgroup }}_j = w_j - \omega \, \sigma_j$, with $\sigma_j$ the standard deviation of $w_j$ and $\omega > 0$ a tuning parameter.
wager2018estimation establish the consistency and asymptotic normality of random forest predictions. However, besides some regularity and technical assumptions, there are certain conditions regarding the construction of individual trees that must be satisfied. In the following, I discuss these conditions.
The first condition requires that the trees use different observations to place the splits and compute the leaf predictions. This condition is called honesty and is crucial to bounding the bias of forest predictions.
wager2018estimation implement honesty by drawing a subsample $\mathcal{S}_b$ from the original sample $\mathcal{S}$ and splitting the subsample into two halves $\mathcal{S}_b^{tr}$ and $\mathcal{S}_b^{hon}$, using $\mathcal{S}_b^{tr}$ to grow the $b$-th tree and $\mathcal{S}_b^{hon}$ to compute its leaf predictions athey2019generalized. Alternatively, lechner2022modified propose a different approach. They divide the original sample $\mathcal{S}$ into a training sample $\mathcal{S}^{tr}$ and an honest sample $\mathcal{S}^{hon}$, constructing trees from random subsamples of $\mathcal{S}^{tr}$ and computing their leaf predictions from $\mathcal{S}^{hon}$. This strategy ensures that, under i.i.d. sampling, the weights assigned to individual units in $\mathcal{S}^{hon}$ are independent of the outcomes of other units, thus allowing for weight-based inference about leaf predictions and their transformations, such as marginal effects. OCF adopts this strategy as well (details in Section (ref)). However, this strategy is somewhat less efficient than the approach proposed by wager2018estimation. This is because, under the latter approach, each data point $w$ will participate in both $\mathcal{S}_b^{tr}$ and $\mathcal{S}_b^{hon}$ of some trees, thus achieving honesty while making more efficient use of the data. However, under this approach, each weight can depend on other units' outcomes, which limits the usage of the weight-based representation of random forest predictions for obtaining standard errors for the leaf predictions and their transformations.
The second condition is that the leaves of the trees must become small in all dimensions of the covariate space as the sample size increases. This is necessary for achieving consistency of the predictions and is accomplished by introducing randomness in the tree-growing process and enforcing a regularity condition on how quickly the leaves get small.
To achieve $\alpha$-regularity, OCF ignores splits that do not satisfy this condition. The algorithm always selects the best split from among the candidate splits that would maintain at least a fraction $\alpha$ of the parent node's observations on both sides of the split. This way, we can rule out any influence of the splitting rule on the shape of the final leaves.
Third, trees must be constructed using subsamples drawn without replacement, rather than bootstrap samples, as originally proposed by breiman2001random.
Fourth, to establish asymptotic normality, trees must be symmetric.
Under these conditions, wager2018estimation establish consistency and asymptotic normality of the random forest predictions. If the $M$ forests constructed by OCF satisfy these conditions, then they inherit these properties, thus producing consistent and asymptotically normally distributed predictions of conditional probabilities.
In addition to the consistency and asymptotic normality of the random forest predictions, wager2018estimation show that the asymptotic variance of such predictions can be consistently estimated by adapting the infinitesimal jackknife estimator proposed by wager2014confidence to the case of subsampling without replacement. This approach can be used to estimate the variance of a prediction $\hat{p}_m^{OCF} \left( \cdot \right)$ at $w$. However, generalizing this method to estimate the variance of marginal effects $\nabla^j \hat{p}_m^{OCF} \left( \cdot \right)$ is not straightforward.
To overcome this limitation, OCF employs an alternative approach that leverages the weight-based representation of random forest predictions athey2019generalized and adapts the weight-based inference proposed by lechner2022modified lechner2019random. In particular, OCF implements honesty in a way that guarantees that the weight assigned to the $i$-th unit is independent of the outcomes of other units. This allows for the derivation of a straightforward formula for the variance of honest predicted probabilities and marginal effects.
First, we express OCF predictions as weighted averages of the outcomes. Let $\mathcal{S}$ denote the observed sample. The following provides an expression for a prediction $\hat{p}_m^{OCF} \left( \cdot \right)$ at $w$ numerically equivalent to that in ((ref)):
where the weights $\hat{\alpha}_{m, 1} \left( w \right), \dots, \hat{\alpha}_{m, |\mathcal{S}|} \left( w \right)$ determine the forest-based adaptive neighborhood of $w$. They represent how often the $i$-th observation in $\mathcal{S}$ shares a leaf with $w$ in the $m$-th forest. This measures how important the $i$-th observation is for fitting $p_m \left( \cdot \right)$ at $w$. Notice that $\sum_{i \in \mathcal{S}} \hat{\alpha}_{m, i} \left( w \right) = 1$ for all $w$.
Calculating the variance of a prediction $\hat{p}_m^{OCF} \left( w \right)$ in ((ref)) is challenging because the weight assigned to the $i$-th unit $\hat{\alpha}_{m, i} \left( w \right)$ is a function of both $\mathcal{S}$ and $W_i$. Thus, this weight depends on the outcomes of all other units in $\mathcal{S}$, which complicates the formula for the variance.
However, the formula for the variance simplifies under the particular honesty implementation of OCF. Let $\mathcal{S}^{tr}$ and $\mathcal{S}^{hon}$ be a training sample and an honest sample obtained by randomly splitting the observed sample $\mathcal{S}$. Also, let $\hat{\alpha}_{m, i}^{tr} \left( \cdot \right)$ be the weights induced by a forest constructed using only $\mathcal{S}^{tr}$. Then, an honest prediction $\tilde{p}_m^{OCF} \left( \cdot \right)$ at $w$ is obtained by the following weighted average of observations in $\mathcal{S}^{hon}$:
The new weight assigned to the $i$-th unit $\hat{\alpha}_{m, i}^{tr} \left( w \right)$ is a function of $\mathcal{S}^{tr}$ and of $W_i$. Thus, under i.i.d. sampling this weight is independent of the outcomes of other units in $\mathcal{S}^{hon}$. This allows us to derive a simple formula for the variance of an honest prediction $\tilde{p}_m^{OCF} \left( w \right)$:
We can estimate this variance by its sample analog.
By plugging ((ref)) into ((ref)), we obtain the following estimator of honest marginal effects:\footnote{\ Similar results are obtained for discrete covariates by plugging ((ref)) into ((ref)).}
with $\check{\alpha}_{m, i}^{tr} ( \widehat{\lceil w_j \rceil}, \widehat{\lfloor w_j \rfloor} ) = \hat{\alpha}_{m,i}^{tr} ( \widehat{\lceil w_j \rceil} ) - \hat{\alpha}_{m, i}^{tr} ( \widehat{\lfloor w_j \rfloor} )$ a transformation of the original weights. Using the same argument as before, under i.i.d. sampling the weight assigned to the $i$-th unit $\check{\alpha}_{m, i}^{tr} ( \widehat{\lceil w_j \rceil}, \widehat{\lfloor w_j \rfloor} )$ is independent of the outcomes of other units in $\mathcal{S}^{hon}$. Thus the variance of an honest marginal effect $\nabla^j \tilde{p}_m^{OCF} \left( w \right)$ can be expressed as follows:
As before, we can estimate this variance by its sample analog.
Following the discussion of Section (ref), the honest predicted probabilities in ((ref)) are consistent and asymptotically normal, provided that the weights $\hat{\alpha}_{m, i}^{tr} \left( \cdot \right)$ are induced by a forest composed of $\alpha$-regular with $\alpha \leq 0.2$ and symmetric random-split trees grown using subsampling without replacement. With these conditions met, we can use the estimated standard errors of honest predicted probabilities $\tilde{p}_m^{OCF} \left( \cdot \right)$ to conduct valid inference as usual, e.g., by constructing conventional confidence intervals.
Furthermore, under the same conditions the honest marginal effects in ((ref)) are a linear combination of normally distributed predictions, and thus have a normal distribution as well. Therefore, we can also construct conventional confidence intervals for honest marginal effects $\nabla^j \tilde{p}_m^{OCF} \left( \cdot \right)$ using their estimated standard errors.
This section uses synthetic data to evaluate the performance of the ordered correlation forest (OCF) estimator. In the next subsection, I present the DGPs employed in the simulation. Then, I compare OCF with various alternative methods in terms of estimating conditional choice probabilities. Finally, I assess the ability of OCF in estimating and making inference about the covariates’ marginal effects.
Latent outcomes are generated as in ((ref)), with $\epsilon_i \sim \textit{logistic} \left( 0, 1 \right)$. Six raw covariates are generated as $W_{i, 1}, W_{i, 3}, W_{i, 5} \sim \mathcal{N} \left( 0, 1 \right)$ and $W_{i, 2}, W_{i, 4}, W_{i, 6} \sim \textit{Bernoulli} \left( 0.4 \right)$. Covariates are independent of one another and of $\epsilon_i$. We consider $W_{i, 5}$ and $W_{i, 6}$ as “ noise" covariates, as they enter the DPGs below with null coefficients.
I consider three designs that differ in the regression function $g \left( \cdot \right)$:
with $\beta = \left( 1, 1, 1/2, 1/2, 0, 0 \right)$ in all designs. Design 1 represents a linear model where all the raw covariates enter without transformation, serving as a benchmark for assessing the performance of the estimators under a straightforward and interpretable setting. In Design 2, the covariates are transformed while preserving the additive structure of the model, thus allowing us to evaluate the estimators' ability to handle non-linearities arising from covariate transformations. Design 3 introduces more complex non-linearities by departing from the additive model structure and employing a nonlinear regression model. For each design, I consider four sample sizes, $| \mathcal{S} | \in \left\{ 500, \, 1000, 2000, 4000 \right\}$. Thus, I consider overall twelve different scenarios.
In each design, I obtain the observed outcomes $Y_i$ by discretizing $Y_i^*$ into three classes:
I construct the threshold parameters $\zeta_1$ and $\zeta_2$ as follows. First, I fix two values $\zeta_1^q = 0.33$ and $\zeta_2^q = 0.66$. Then, I generate a sample of $1,000,000$ $Y_i^*$ and set $\zeta_m = Q \left( \zeta_m^q \right)$, with $Q \left( \cdot \right)$ the empirical quantile function of $Y_i^*$. This way, the threshold parameters are uniformly spaced, and the class widths are approximately equal.
After drawing a sample $\mathcal{S}$, I estimate the conditional choice probabilities using both multinomial and ordered machine learning techniques, combining them with random forests breiman2001random and penalized logistic regressions with an L1 penalty tibshirani1996regression. I refer to the resulting estimators as multinomial random forest ($MRF$), multinomial L1 regression ($ML1$), ordered random forest ($ORF$), and ordered L1 regression ($OL1$). I also consider two versions of OCF, the “ adaptive" version $OCF_{\scriptstyle_{A}}$ and the “ honest" version $OCF_{\scriptstyle_{H}}$. This way, we can quantify the loss in the precision derived from using fewer observations to build the forests, representing the price to pay for valid inference. Finally, I include the standard ordered logit ($LOGIT$) model as a parametric benchmark for comparison.
To account for non-linearities in $g \left( \cdot \right)$, the parametric methods $LOGIT$, $ML1$, and $OL1$ employ different linear-in-parameter models such as ((ref)). Three different specifications are considered. The first specification consists of a model with only the raw covariates, that is, with $X_i = W_i$. The second specification introduces third-order polynomials for continuous covariates, leading to a set of $12$ covariates. The third specification enlarges this set by adding all the two-way interactions between the raw covariates, resulting in a total of $27$ covariates. In contrast, the forest-based estimators $MRF$, $ORF$, $OCF_{\scriptstyle_{A}}$, and $OCF_{\scriptstyle_{H}}$ are fed with only the raw covariates without adding any polynomials, interaction terms, or other transformations of the covariates, as these estimators can naturally handle non-linearities in $g \left( \cdot \right)$. To implement $OCF_{\scriptstyle_{H}}$, I randomly split $\mathcal{S}$ into a training sample $\mathcal{S}^{tr}$ used to construct the trees and an honest sample $\mathcal{S}^{hon}$ used to compute the leaf predictions. I choose $| \mathcal{S}^{tr} | = | \mathcal{S}^{hon} | = | \mathcal{S} |/2$.
I rely on an external validation sample $\mathcal{S}^{val}$ of size $| \mathcal{S}^{val} | = 10,000$ to assess the predictive performance of the estimators. This large number of observations helps minimize the sampling variance. For each replication $r = 1, \dots, R$, I calculate the mean squared error, mean absolute error, and ranked probability score for each estimator:
with $\hat{p}_{m, r} \left( \cdot \right)$ the estimated conditional probabilities in the $r$-th replication, and $\hat{\mu}_{m, r} \left( w \right) = \sum_{j = 1}^m \hat{p}_{j, r} \left( w \right)$ the estimated cumulative distribution function. Notice that, by simulation design, we can compute the true probabilities as in ((ref)). I summarize these performance measures by averaging over the replications.\footnote{\ The objective of this subsection is to evaluate the prediction accuracy of each estimator. Thus, I do not consider the variance or the actual coverage rates of confidence intervals as performance measures, as these aspects are not relevant when the interest lies in prediction accuracy.}
\begingroup {8pt}
\endgroup
Table (ref) displays the results obtained with $R = 2,000$ replications. The simulation shows that OCF outperforms all other forest-based estimators uniformly across all considered scenarios. In particular, OCF consistently achieves lower MSE and MAE than $MRF$ and $ORF$, and minimal disparities in RPS are observed. An exception arises in Design 3 where, for the two smallest sample sizes, OCF and $MRF$ show similar performance.
The simulation also shows that OCF maintains a competitive performance when compared to the parametric estimators $LOGIT$, $ML1$, and $OL1$. As expected, when fed with only the raw covariates, $LOGIT$ and $OL1$ perform best in Design 1 since they correctly specify the parametric model and the distributional assumption of the error term. These estimators are closely followed by $ML1$ in terms of performance. However, the performance of $LOGIT$, $ML1$, and $OL1$ fed with only the raw covariates deteriorates when non-linearities in $g \left( \cdot \right)$ are introduced, causing them to rank among the worst estimators. For the smallest sample size, the performance gap with respect to OCF is relatively substantial in Design 2 (between $26$--$36\%$ in terms of MSE, between $8$--$16\%$ in terms of MAE, and between $43$--$48\%$ in terms of RPS) and moderate in Design 3 (around $11\%$ in terms of MSE and around $17\%$ in terms of RPS, with minimal disparities in MAE observed). However, in larger samples, this performance gap becomes more pronounced, with the MSE of $LOGIT$, $ML1$, and $OL1$ reaching values up to $165\%$ larger than that of OCF, MAE up to $61\%$, and RPS up to $232\%$.
Adding constructed covariates when the model for $g \left( \cdot \right)$ is linear and the raw covariates enter without transformation (Design 1) deteriorates the performance of $LOGIT$ and $OL1$, as this primarily inflates the variance of the estimation. This effect becomes less relevant in larger samples. In contrast, when non-linearities in $g \left( \cdot \right)$ are introduced - either via transformations of the covariates (Design 2) or by employing a non-linear regression model (Design 3) - adding polynomials and interactions of the covariates significantly improves the performance of the parametric estimators. In particular, in Design 2, $LOGIT$, $ML1$, and $OL1$ achieve their best performance by introducing third-order polynomials, with their performance deteriorating when interactions between covariates are also included, although this deterioration is substantially attenuated in larger samples. However, in Design 3, where more complex non-linearities are introduced, including the interactions is necessary to achieve the best possible performance.
When constructed covariates are included in their specifications, $LOGIT$, $ML1$, and $OL1$ generally exhibit lower MSE, MAE, and RPS than OCF. However, in Design 2, OCF outperforms all the parametric methods in larger samples, with advantages over the best parametric specification ranging between $22$--$31\%$ in terms of MSE, $7$--$10\%$ in terms of MAE, and $47$--$55\%$ in terms of RPS. An exception arises in the largest sample where $LOGIT$ and OCF tie in terms of MAE. Moreover, in Design 3, the performance gap appears to narrow as the sample size increases, suggesting that OCF might outperform the parametric methods if enough observations are used to train the model.
Finally, we compare the adaptive and the honest versions of OCF to quantify the price to pay for valid inference. Surprisingly, the honest version $OCF_{\scriptstyle_{H}}$ performs better than the adaptive version $OCF_{\scriptstyle_{A}}$ in almost all scenarios despite using half of the observations to construct the forests. Honesty reduces the bias of the forests' estimates but generally comes at the expense of a higher variance. In this simulation, the reduction in bias appears to outweigh the increase in variance, resulting in improved prediction performance.
After drawing a sample $\mathcal{S}$, I split it into a training sample $\mathcal{S}^{tr}$ and an honest sample $\mathcal{S}^{hon}$ of equal size. Then, I use $\mathcal{S}^{tr}$ to construct the forests, and $\mathcal{S}^{hon}$ to estimate honest marginal effects at the mean and median of the covariates as in equation ((ref)). Additionally, I use the sample analog of equation ((ref)) to get standard errors for the estimated effects.\footnote{\ Estimating the mean or the median marginal effect and its standard error would involve computing the weights $\check{\alpha}_{m, i}^{tr} \left( \cdot, \cdot \right)$ for each prediction point $w$, which would result in an impractically long computational time for a Monte Carlo exercise. Therefore, I restrict the analysis solely to the marginal effects at the mean and median of the covariates.}
To assess the performance of the estimator, I calculate the squared bias and variance for each marginal effect, as well as the actual coverage rates of their corresponding $95\%$ confidence intervals. Notice that, by simulation design, we can compute the true marginal effects. I summarize these performance measures by averaging across all marginal effects.
Table (ref) displays the results obtained with $2,000$ replications. Overall, the simulation shows the ability of OCF to conduct asymptotically valid inference about marginal effects. The estimated squared bias consistently remains close to zero, indicating that the estimator is approximately unbiased. In smaller samples, the actual coverage rates of the confidence intervals tend to fall below the nominal rate and can be as low as $70\%$. However, as the sample size increases, the coverage rates gradually converge to the nominal level. In Design 2 and Design 3, more observations are required to reach the nominal rate compared to Design 1.
\begingroup {8pt}
\endgroup
When we compare these results with those presented in Table (ref), an interesting pattern emerges: the actual coverage rates of the confidence intervals tend to be worse when the predictive performance of $OCF_{\scriptstyle_{H}}$ is lower. This pattern aligns with the fact that OCF estimates marginal effects by post-processing its conditional probability predictions (see equations (ref)--(ref)), meaning that the quality of conditional probability estimation directly impacts the accuracy of marginal effects estimation. As evident in Table (ref), $OCF_{\scriptstyle_{H}}$ performs best in Design 1 and exhibits relatively lower performance in Design 3 compared to Design 2. Consequently, for any given sample size, we observe a larger estimated squared bias in Design 2 and Design 3 relative to Design 1, which explains why more observations are needed to reach the nominal rate in these designs. However, as the sample size increases, the predictive performance of OCF improves, and thus the estimated bias decays asymptotically. Therefore, in larger samples, the estimated confidence intervals are more likely to be centered around the true estimand, resulting in actual coverage rates that converge to the nominal level.
This section uses real data to compare the predictive performance of the ordered correlation forest estimator with the same estimators of Section (ref).
\begingroup {8pt}
\endgroup
I utilize the same data sets considered by janitza2016random, hornung2020ordinal, and lechner2019random. These data sets differ in terms of the number of covariates, observations, and classes of the observed outcome. Table (ref) provides a summary of the data sets. For further details on the background of each data set, the reader is referred to janitza2016random.
To assess the prediction accuracy of each estimator, I employ a ten-fold cross-validation procedure. Specifically, I randomly divide each data set into ten folds $\mathcal{S}^1, \dots, \mathcal{S}^{10}$ with roughly equal sizes. For each fold $f = 1, \dots, 10$, I fit all the estimators using the observations from all the other folds except for $\mathcal{S}^f$. Then, I calculate the same performance measures of Section (ref) using the held-out $\mathcal{S}^f$:
with $\hat{p}_{m, f} \left( \cdot \right)$ the estimated conditional probabilities using all the other folds except for $\mathcal{S}^f$, and $\hat{\mu}_{m, f} \left( w \right) = \sum_{j = 1}^m \hat{p}_{j, f} \left( w \right)$ the estimated cumulative distribution function. Finally, I repeat this process ten times. This approach eliminates the dependence of the results on a particular training-validation sample split.
Figure (ref) reports the results by displaying boxplots showing the median and interquartile range of the estimated MSE, MAE, and RPS, together with their minima and maxima.\footnote{\ The cross-validation exercise yields a smaller sample size compared to the simulation results presented in Section (ref). Consequently, estimates of expected MSE, MAE, and RPS can be more imprecise and influenced by outliers. I report the distribution of the estimated MSE, MAE, and RPS using boxplots to provide a more robust assessment of the prediction performance of each estimator.}
Overall, the results indicate that OCF performs competitively compared to the other estimators, with no substantial differences in performance observed in most data sets. In the smallest data set under consideration (vlbw), $LOGIT$, $ML1$, and $OL1$ perform marginally better than $MRF$, $ORF$, and OCF in terms of MSE, with similar MAE and RPS observed. However, in the largest data set (wines), OCF emerges as one of the best estimators together with $MRF$ and $ORF$. This result highlights the advantage of forest-based methods over parametric methods in larger samples.
The addition of constructed covariates deteriorates the performance of $LOGIT$ while it does not substantially change the performance of $ML1$ and $OL1$, except in the largest data set (wines) where it leads to improved performance for all parametric estimators. The deterioration in the performance of $LOGIT$ is particularly pronounced in the support and nhanes data sets. In these data sets, including third-order polynomials for continuous covariates and all the two-way interactions between the raw covariates results in a total of $324$ and $1394$ covariates, respectively. Given the moderate sample sizes, the inclusion of the additional covariates substantially increases the variance of the estimation, causing $LOGIT$ to perform worse compared to a specification with only the raw covariates. In contrast, the performance of $ML1$ and $OL1$ is not substantially affected by the inclusion of the additional covariates, as these estimators employ regularization techniques that help contain the increase in variance.
In contrast to the simulation results, the honest version of OCF does not outperform the adaptive version. In the two smaller data sets (vlbw and mammography), $OCF_{\scriptstyle_{A}}$ and $OCF_{\scriptstyle_{H}}$ exhibit similar MSE and RPS, with the MAE of $OCF_{\scriptstyle_{H}}$ being slightly larger than that of $OCF_{\scriptstyle_{A}}$. However, as the sample size increases, a performance gap in favor of $OCF_{\scriptstyle_{A}}$ emerges in terms of all performance measures. This gap becomes substantial in the largest data set.
This paper proposes a novel machine learning estimator specifically optimized for handling ordered non-numeric outcomes. The proposed estimator adapts a standard random forest splitting criterion breiman2001random to the mean squared error relevant to the specific estimation problem at hand, thus mitigating the biases that traditional methods can introduce. The new splitting rule is then used to build a collection of forests, each estimating the conditional probability of a single class. A nonparametric approximation of derivatives is employed to estimate the covariates' marginal effects lechner2019random.
Under an “ honesty" condition athey2016recursive, the estimator inherits the asymptotic properties of random forests, namely the consistency and asymptotic normality of their predictions wager2018estimation. The particular honesty implementation used by the ordered correlation forest allows us to obtain standard errors for the covariates' marginal effects by leveraging the weight-based representation of the random forest predictions athey2019generalized. The estimated standard errors can then be used to construct asymptotically valid symmetric confidence intervals.
Evidence from synthetic data shows that the proposed estimator features a superior prediction performance than alternative forest-based estimators and demonstrates its ability to construct valid confidence intervals for the covariates' marginal effects.
\nocite{*}
\singlespacing \newrefcontext[sorting = nty] \printbibliography