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.
75,560 characters · 15 sections · 46 citation commands
Asymptotic Normality for Multivariate Random Forest Estimators
\onehalfspacing
Trees and random forests are non-parametric estimators first introduced by Breiman breiman01. Given a feature space $\mathcal{X} \subset \mathbf{R}^p$ and a set of data points $\{ (X_i, Y_i) \} \subseteq \mathcal{X} \times \mathbf{R}$, tree estimators recursively partitions the feature space into axis-aligned non-overlapping hyperrectangles\footnote{When the feature space $\mathcal{X}$ need not be rectangular, one may always enlarge $\mathcal{X}$ to a rectangular set $\mathcal{X}'$ that is defined to the intersection of all rectangular sets containing $\mathcal{X}$.} by repeatedly splitting $\mathcal{X}$ along a given axis. The prediction of the tree estimator at a test point $x \in \mathcal{X}$ is then an aggregate of the targets $Y_i$'s that land in hyperrectangle containing $x$; when $Y_i$ is continuous, the aggregate is the sample mean and the tree is also known as a regression tree. The depth of a tree estimator---defined as the maximal number of splits taken before reaching a terminal hyperrectangle---controls the complexity of the tree estimator. There are two popular methods for controlling complexity: the “boosting” approach grows trees of large depth, then reduces complexity by either trimming the tree (i.e., so that predictions are made at a non-terminal hyperrectangle) or introducing a decay factor; the “bagging” approach instead grows a collection of shallow trees on different subsets of the data, and averages over the trees for the final prediction. The intuition for bagging is that trees grown on different subsets are not perfectly correlated, so that aggregation reduces variance and balances the bias-variance tradeoff. Estimators of this type are called random forests, and they are the focus of this paper.
Since their introduction in the early 2000s, random forests have become an increasingly important tool in applied data analysis, owing to a multiple of practical advantages over competing models. First, high-quality random forest libraries are readily available, with popular implementations that scale to hundreds of distributed workers lightgbm,xgboost. Moreover, the core algorithm behind tree estimators and random forests are simple enough to allow for rapid prototyping of bespoke implementations, e.g. grfgithub. Another advantage of tree-based methods is that they can ingest real-world data without much issue: continuous, discrete, and ordered categorical features may be freely mixed\footnote{Splits on discrete features partitions that variable into two arbitrary non-empty sets; no changes are needed for ordered categorical features.} catboost, model estimates are immune to feature outliers, and missing data may be easily incorporated. Firstly, their construction naturally aligns with the spatial locality of most applied: that is, the underlying target function relating $Y$ to $X$ is continuous. Finally, tree models are interpretable, with well-defined notions of feature importance gregorutti2017correlation,strobl2008conditional, which supports their use as model selection tools genuer2010variable.
Within economics, random forests may be fruitfully applied to estimate heterogeneneous treatment effects. In Rubin's potential outcomes framework rubin (see imbens_rubin_2015 for an overview), an individual $i$ is associated with two potential outcomes $Y^{(0)}$ and $Y^{(1)}$, with one of the outcomes being realized depending on whether $i$ undergoes treatment. The statistician has access the IID observations $\{ (X_i, W_i, Y_i) : 1 \leq i \leq n \}$, where $X_i$ is a vector of observed covariates for individual $i$, $W_i \in \{ 0, 1 \}$ is (an encoding of) their treatment status, and $Y_i = Y_i^{(W_i)}$ is her realized outcome. One quantity of interest is the treatment effect at $x$
Since only one of $Y^{(0)}_i$ and $Y^{(1)}_i$ is observed, consistent estimation of $\tau(x)$ requires further distributional assumptions. A common assumption is unconfoundedness, i.e., that treatment status $W_i$ is independent of $Y^{(1)}$ and $Y^{(0)}$ conditional on $X_i$. Under this assumption,
Here, the key function is $e(x)$, known as the propensity score, is the probability of treatment for the subpopulation with covariates $x$; see propensity_score the derivation and implications. Machine learning methods---including random forests---may be brought to bear on the problem by estimating $e(x)$. Alternatively, unconfoundedness also implies
so that $\tau(x)$ may be estimated by fitting two models, one on the subset of the sample in which $W = 1$, and the other on $W = 0$.
In econometric applications, conducting pointwise inference on the target function $f: \mathcal{X} \to \mathbf{R}$ (e.g., to test the null hypothesis $H_0: f(x) = 0$) requires knowledge about about the rate of convergence or asymptotic distribution of the underlying estimator $\hat f(x)$, where $x$ is the point of interest. However, functionals of target function are often also of interest: for example, the difference of treatments effects (i.e., $f = \tau$) for two different subpopulations is captured by the quantity
where $x$ and $\bar x$ are covariates describing the two subpopulations. More generally, we might also be interested in a weighed treatment effect, where a subpopulation $x$ is given an importance weight modeled as a density $\mu(x)$. In this case, the corresponding functional of $f$ is
and the integral is taken over the domain $\mathcal{X}$.
Inference on functionals of $f$ requires not only the asymptotic distribution of the point estimate $f(x)$, but also the correlation between estimates at different points $f(x)$ and $f(\bar x)$. As a concrete example, consider the function $\tau(x)$ and the simple difference $\tau(x) - \tau(\bar x)$. We have
We may estimate the difference by estimating $A$ and $B$ separately, fitting a random forest model to the two “halves” of the dataset where $W_i = 1$ and $W_i = 0$, as discussed above. The estimators $\hat A$ and $\hat B$ obtained are thus independent, so that $\operatorname{\mathbf{Var}}(\hat A - \hat B)=\operatorname{\mathbf{Var}} \hat A + \operatorname{\mathbf{Var}} \hat B$. The variances $\hat A$ and $\hat B$ then depend on the covariance of their respective random forest estimates at $x$ and $\bar x$.
This paper studies the correlation structure of a class of random forests models whose asymptotic distributions were first worked out in crf. We find sufficient conditions under which the asymptotic covariance of random forest estimates at different points vanish relative to their respective variances; moreover, we provide finite sample heuristics based on our calculations. To the best of our knowledge, this is the first set of results on the correlation structure of random forest estimators.
The present paper builds on and extends the results in crf, which in turn builds on related work Wager2015AdaptiveCO on general concentration properties of trees and random forest estimators. See also athey2019, which extends the random forest model considered here to a broader class of target functions by incorporating knowledge of moment conditions. Stability results established in this paper have appeared in arsov2019stability, who study the notions of algorithmic stability for random forests and logistic regression and derive generalization error guarantees. Also closely related to our paper are chernozhukov2017 and chen2018, concerning finite sample Gaussian approximations of sums and $U$-statistics in high dimensions. In this context, our paper provides a stepping stone towards applying the theory of finite sample $U$-statistics to random forests, where bounds on covariance matrices plays a central role.
The paper is structured as follows. In Section 2, we introduce the random forest model and state the assumptions required for our results; Section 3 contains our main theoretical contributions; Section 4 builds on Sections 3 and discusses heuristics useful in finite sample settings; Section 5 concludes. All proofs are found in the appendix.
The goal of this paper is to study asymptotic Gaussian approximations of random forest estimators. Throughout, we assume that a random sample $\{ Z_i = (X_i, Y_i) : 1 \leq i \leq n \} \subseteq \mathcal{X} \times \mathbf{R}$ is given, where each $X_i$ is a vector of features or covariates belonging to a subset $\mathcal{X} \subseteq \mathbf{R}^p$ of $p$-dimensional Euclidean space, and $Y_i \in \mathbf{R}$ is the response or target corresponding to $X_i$. We will refer to $\mathcal{X}$ as the feature space or the feature domain.
Given the data set $\{ Z_i \}_{i=1}^n$, a tree estimator recursively partitions the feature space by making axis aligned splits. Specifically, an axis-aligned split is a pair $(j, t)$ where $j \in \{ 1, \dots, p \}$ is the splitting coordinate and $t \in \mathbf{R}$ is the splitting index; given a subset $\mathcal{R} \subseteq \mathcal{X}$, a split $(j, t)$ divides $\mathcal{R}$ into left and right halves
where $x_j$ denotes the $j$-th coordinate of the vector $x$. Starting with the entire feature space $\mathcal{X}$, the recursive splitting algorithm computes a (axis-aligned) split based on the data $\{ Z_i: 1 \leq i \leq n \}$; for example, when the target $Y_i$ continuous, a popular choice is
where $L = L(\tilde j, \tilde t)$ and $R = R(\tilde j, \tilde t)$ are the two halves of $\mathcal{X}$ obtained by the split $(\tilde j, \tilde t)$, with $\mu_L$ and $\mu_R$ being the averages of targets $Y_i$ whose corresponding feature land in $L$ and $R$, respectively.
After the first split, $\mathcal{X}$ is split into two halves $L$ and $R$. The process is then repeated for $L$ and $R$ separately, in that a split for $L$ is computed by using the subset of the data whose features $X_i$ belong in $L$, and likewise for $R$. Each of the halves is then split again, and so on, until a stopping criterion is met. The process completes when the stopping criterion is satisfied for each node; at this point, the collection of halves form a partition\footnote{According to (ref), we exclude edge cases where points lead on the an “edge” of the rectangle. This is not an issue for continuous variables, while for categorical variables, the definition should be slightly changed so that one of the halves contain $x_j = t$. In this paper, we deal only with continuous features.} of $\mathcal{X}$, with each partition---and all the halves that came before it---being the intersection of a hyperrectangle with $\mathcal{X}$. The sequence of splits corresponds to a tree in the natural way; we will call the halfspaces that arise during the splitting as nodes, and elements of the final partition terminal nodes.
Given the collection $N_1, \dots, N_q$ of terminal nodes (which form a partition of $\mathcal{X}$), the prediction of the tree at generic test point $x \in \mathcal{X}$ is the average of the responses that belong in the same terminal node as $x$
where the outer sum runs over observations $i$ for which $X_i$ belongs to the partition $N_j$, and $|N_j|$ is the number of such observations. The input $\xi$ is an external source of randomization to allow for randomized split selection procedures. Thus, $T(x; \xi, Z_1,\dots,Z_n)$ refers to prediction at $x$ for a tree grown using data $\{ Z_1, \dots, Z_n \}$ with randomization parameter $\xi$. As a function of $x$, keeping $Z_1, \dots, Z_n$ and $\xi$ fixed, $T:\mathcal{X} \to \mathbf{R}$ is then a step function, i.e., a linear combination of indicator functions of rectangular sets.
We note here that equations (ref) and (ref) are not the only possible choices; in particular, the rule used to choose the optimal split may be path dependent (i.e., dependent on previous splits) as in popular implementations (see xgboost, which allows for random forests but uses gradient boosting after the initial split), and the final prediction rule (ref) may instead do a final linear fit or use weighted average (c.f., lightgbm,xgboost). Analysis of tree (and random forest) models is complicated in these cases, so the present paper stipulates that the algorithm uses a splitting rule that is “similar” to similar to (ref) (c.f. Proposition (ref)) and uses (ref) as the final prediction. In particular, this implies that the target estimated by the tree estimator is the regression function $x \mapsto \operatorname{\mathbf{E}}(Y \mid X = x)$.
Given a specific tree estimator $T$, i.e., given a set of splitting rules which defines the tree estimator, we define the random forest estimator to be the average of the tree estimator across all $\binom{n}{s}$ subsamples, marginalizing over the randomization device $\xi$. Specifically, the random forest estimate $\operatorname{RF}(x)$ at $x \in \mathcal{X}$, given data $\{ Z_1, \dots, Z_n \}$, is defined to be
where the summation runs over size-$s$ subsets of $\{ 1, \dots, n \}$, and the inner expectation is taken with respect to $\xi$. Importantly, note each tree is grown on a subsample of size $s<n$. We follow crf in assuming that $s \sim n^{\beta}$ for some $\beta$ sufficiently close to one; specifically, we assume throughout that the subsample size is chosen as to satisfy the assumptions of Theorem 3 of crf, so that---along with other assumptions to be introduced presently---the random forest estimator $\operatorname{RF}$ is a consistent estimator of the target function $x \mapsto \operatorname{\mathbf{E}}(Y \mid X = x)$ (see discussion above).
In keeping with the notation crf, we will write $T(x; Z_1, \dots, Z_s)$ to refer to the expectation of $T(x; \xi, Z_1, \dots, Z_s)$ over $\xi$. With this notation, the random forest (at $x$) estimator (ref) is a $U$-statistic with the size $s$ kernel $(Z_1, \dots, Z_s) \mapsto T(x, Z_1, \dots, Z_s)$. We will discuss the $U$-statistic representation of $\operatorname{RF}$ in greater detail in Section 3.
As our results will be an extension of the results in crf, we will study the same model of random forests and adopt a similar set of assumptions. The assumptions regarding tree estimators have appeared before in Wager2015AdaptiveCO, while the distributional assumptions on the conditional moments of $Y$ are standard (see e.g., Chapters 7 and 9 in hastie2009elements).
The first---and most bespoke---assumption is that the tree algorithm is honest. Intuitively, honesty stipulates that that knowledge of the tree structure does not affect the conditional distribution of tree estimates when the features are fixed.
There are several ways to satisfy this assumption. The first is to calculate splits based only the features $X_i$ only. This rules out out the example splitting rule given in (ref), so we could instead use its analog in feature space
where here $\mu_L$ and $\mu_R$ denote the average (i.e., center of mass) of the $X_i$ in each halfspace. In this instance, the choice of splits is essentially a clustering algorithm that finds the best division of the sample points into two parts. Another way to satisfy to the honesty assumption while still computing splits based on the targets is to use sample splitting. The data is partitioned into two parts $\mathcal{I}_1$ and $\mathcal{I}_2$; observations in $\mathcal{I}_1$ and $X_i \in \mathcal{I}_2$ may be freely used during the splitting process, while $Y_i \in \mathcal{I}_2$ are used to determine terminal node values. In this case, equality in (ref) is required to hold for $i \in \mathcal{I}_2$. Finally, a third method satisfy honesty requires the existence of auxiliary data $\{ W_i \}$. During the splitting stage (“model fitting”), splits are computed as if the response variable is $W_i$; for example, (ref) is used with $\mu_L$ and $\mu_R$ being the average of the $\{ W_i \}$. Once the tree is fully grown, predictions (“model inference\footnote{The terminology `inference' is used to mean computing the predictions of an existing model, which is unrelated to the typical usage of `inference' in econometrics. The former terminology is standard in applied settings, used when describing a data pipeline: see the documentation of tensorflow, sklearn.}”) are made using $Y_i$'s as usual. The practice of using such surrogate targets is especially popular in time-series prediction, where different horizons are used in fitting and inference steps (c.f., quaedvlieg2019multi).
In the present paper, for simplicity of notation, we shall assume that the first scheme is used to satisfy honesty---namely, splitting decisions are based on the feature vectors $X_i$ only. Our results extend to all three schemes.\footnote{ As in crf, constants appearing in our bounds may change in scheme two.}
Our next assumption will ensure that each one of the $p$ axes is chosen as the splitting coordinate with a probability bounded from below.
This is a modification of the random splitting assumption in crf, in which each of the $p$ axes has a probability $\delta$ of being chosen at each split. This could be directly implemented by flipping a $p\delta$ and selecting one of the $p$ coordinates uniformly at random to split when the coin lands heads. Another method, studied in athey2019 and implemented by popular libraries such as xgboost, uses randomizes the number of available splitting axes in each round. Specifically, a Poisson random variable $Q$ with intensity proportional to $\sqrt{p}$ is first realized ($Q$ is realized independently from round to round). Afterwards, $\min(Q, p)$ many axes are uniformly selected as potential candidates\footnote{Splitting for that node ceases if $Q = 0$.} for splitting in that round. Clearly, this also leads to a lower bound on the probability of being chosen for each coordinate $1 \leq j \leq p$.
Each of the two methods above involves two separate rounds of randomization: first, a random variable encoding the decision to split randomly is made (i.e., the coin or the random variable $Q$), and second, the splitting axis is then determined. Intuitively, our cyclic splitting assumption above forgoes the second randomization step: in doing so, the variance of the number of times that any coordinate is chosen is reduced. Importantly, the variance will depend only on $\delta$ and not on $p$, which we will exploit in our proofs.
This key assumption is carried over from crf and contains two requirements. The first requires that no split may produce a halfspace containing too few observations, i.e., that both halfspaces are large when measured the by count of observations. As shown in Wager2015AdaptiveCO, this implies that with exponentially small complementary probability, the splitting axis shrinks by a factor between $\alpha$ and $1-\alpha$, so that both halfspaces also large in Euclidean volume (with high probability).
The second half of the assumption places an upper bound on the number of observations in terminal nodes. Trees grown under this assumption will necessarily be deeper as the sample size $n$ (and thus, the subsample size $s = n^{\beta}$) increases. In particular, the predictions at leaf nodes---averages of observations $Y_i$--- will be averages of a bounded number of terms. An important consequence is that the variance of the tree estimator (at any test point $x$) is bounded below (c.f., the distributional assumption on $\operatorname{\mathbf{Var}}(Y \mid X = x)$).
The predetermined splitting assumption is specific to our paper. That candidate splits considered at each node are data-independent is “almost” without loss of generality since the feature space $\mathcal{X}$ is fixed. For example, implementations of random forests typically use 32-bit floating point numbers as their splitting index, so that the assumption is automatically satisfied. Furthermore, we note that the assumption allows candidates splits to depend on the node, so that the splitting process is still data-driven insofar as the sequence of splits leading up to node depends on the observations.
One interpretation of this assumption is that it aids in data compression. For example, suppose that all features are continuous and $\mathcal{X} = [0,1]^p$, and that all candidate splits have the form $(j, k/2^{m})$ for some integers $1 \leq j \leq p$ and $0 \leq k < 2^m$, where $m \geq 1$ is a fixed integer. If a splitting rule such as (ref) is used, then the optimal split depends only on the values $\{ \lfloor 2^m X_{ij} \rfloor : 1 \leq i \leq n, 1 \leq j \leq p \}$ instead of $\{ X_{ij} \}$. Each member of the former set is an integer in $\{ 0, \dots, 2^m-1 \}$ and thus could be represented using $m$ bits. In particular, for a grid resolution of $2^8 = 256$---a fine grid even in moderately large dimensions---each coordinate of the feature vectors $X_i$ may be stored in a single byte. Since modern CPUs and graphics processors store floating point numbers in four or eight bytes, this allows for a substantial reduction, allowing computation power to scale to to larger datasets. The process of encoding in this way is known as quantizing, which is an option supported by many popular software packages. In this way, though the predetermined split assumption may seem at first glance restrictive, it aligns our model more closely with practice.
The continuity and variance bound assumptions are standard. Note that a consequence of continuity and compactness of the hypercube is that the conditional moments up to order three are bounded. Our results will not explicitly depend on knowledge of the density of $X$: however, the density will affect the implicit constants that we carry throughout our proofs (c.f., Lemma 3.2 and Theorem 3.3 in crf).
We begin investigation of the random forest estimator in this section. As discussed in the model introduction, the random forest estimator $\operatorname{RF}(x)$ at a test point $x$ is a $U$-statistic where the kernel is tree estimator $T(x)$ marginalized over external randomizations. This paper studies the multivariate distribution of $\operatorname{RF}$, specifically the correlation structure between $\operatorname{RF}(x)$ and $\operatorname{RF}(\bar x)$ at distinct points $x$ and $\bar x \in \mathcal{X}$. Towards that end, we shall fix a collection of $q$ test points
throughout the remainder of the paper. As these points will remain fixed, for notational brevity we will omit their explicit dependence when writing estimators. Therefore, $\operatorname{RF}(Z_1, \dots, Z_n)$ stands for the $q$-dimensional estimator that is the random forest evaluated at $x_1, \dots, x_q$, given observations $\{ Z_i : 1 \leq i \leq n \}$. As a consequence of notation, most of our equations are to be understood in $\mathbf{R}^q$, with equality and and arithmetic acting coordinate-wise. Finally, when there is no confusion, subscripts (typically $k$ or $\ell$) typically denote a specific coordinate, i.e., the estimate at the $k$-th or $\ell$-th test point; a notable exception is $T_1$, which refers to a Hajek projection that we now describe.
We start by reviewing properties of the H\"oeffding Decomposition of $U$-statistics, also known as Hajek projections; see vdv for a textbook treatment of the univariate case. Let $f(Z_1, \dots, Z_m) \in \mathbf{R}^q$ be a generic $q$-dimensional statistic based on $m$ observations. The Hajek projection of $f$ is defined to be
That is, it is the coordinate-wise projection of $f$ to the linear space spanned by functions of the form $\{ g(Z_i) : 1 \leq i \leq m \}$. In particular, when $f$ is symmetric in its arguments and $Z_1, \dots, Z_m$ is an IID sequence, we have
where is the function such that $f_1(Z_1) = \operatorname{\mathbf{E}}(f \mid Z_1)$, i.e., $f_1(z) = \operatorname{\mathbf{E}}(f \mid Z_1 = z)$.
In our setting, applying the Hajek projection to the centered statistic $\operatorname{RF} - \mu$, where $\mu$ is the expectation of $\operatorname{RF}$, yields
where $i_1, \dots, i_s$ run through the $\binom{n}{s}$ size-$s$ subsets of $\{ 1, \dots, n \}$ as usual. (Recall that $\operatorname{RF}$, $\mu$, and $T$ are all vectors in $\mathbf{R}^q$.) Since the samples $Z_1, \dots, Z_n$ are independent, $\operatorname{\mathbf{E}}(\operatorname{\mathbf{E}}_{\xi}T(\xi;Z_{i_1}, \dots, Z_{i_s}) \mid Z_i) = \mu$ whenever $i \notin \{ i_1, \dots, i_s \}$. As $\{ i_1, \dots, i_s \}$ runs over the size-$s$ subsets of $\{ 1, \dots, n \}$, there are exactly $\binom{n-1}{s-1}$ many which contain $i$. For each of of these subsets,
where $T_1(z) = \operatorname{\mathbf{E}}_{\xi, Z_2, \dots, Z_s} T(\xi; z, Z_2, \dots, Z_s)$. Therefore,
The sequence of observations $Z_1, \dots, Z_n$ is assumed to IID, and this property is preserved for the sequence $\{ T_1(Z_i) : 1 \leq i \leq n \}$ of projections. It is easily verified that $\operatorname{\mathbf{E}}(\mathring \operatorname{RF}) = \mu$, and the point of the previous equation is that it expresses the centered statistic $(\mathring \operatorname{RF} - \mu)$ as an average of centered IID terms, scaled by $s$. This will be our main entry point in establishing asymptotic joint normality.
The standard technique in deriving the asymptotic distribution of a $U$-statistic is to establish a lower bound on the variance of its Hajek projection; this is the approach taken by crf and we follow the approach here. Let $V$ be the variance of $\mathring\operatorname{RF}$; using (ref), we have
where $\mathring T$ is the Hajek projection of the statistic $T$ as in (ref), where $T = \operatorname{\mathbf{E}}_{\xi}T(\xi, Z_1, \dots, Z_s) \in \mathbf{R}^q$.
Since the third moment $\operatorname{\mathbf{E}}(Y^3 \mid X = x)$ is assumed to be bounded, conditions for the Lindberg Central Limit Theorem billingsley2008probability easily follows, and applying the triangular CLT, we have the familiar fact
where $0$ is the zero vector in $\mathbf{R}^q$ and $I$ the $q \times q$ identity matrix.
The asymptotic normality of the random forest estimator $\operatorname{RF}$ can be related to the asymptotic normality of $\mathring \operatorname{RF}$ via
Since the second summand on the RHS is asymptotically normal, by Slutsky's Theorem, $V^{-1/2}(RF-\mu)$ is asymptotically normal once we establish the convergence $V^{-1/2}(\operatorname{RF} - \mathring \operatorname{RF}) \xrightarrow{\;\operatorname{\mathbf{P}}\;} 0$. The strategy is to show that $e = V^{-1/2}(\operatorname{RF} - \mathring \operatorname{RF})$ converges in squared mean. We may develop its squared norm via
where we used the identity $\operatorname{tr}(ABC) = \operatorname{tr}(BCA)$ for conforming matrices $A$, $B$, and $C$. That trace on the extreme RHS goes to zero is the natural multivariate generalization of the familiar condition
for univariate $U$-statistics vdv.
In the univariate setting, the previous condition is checked by considering higher order decompositions proceeds by considering higher order decompositions of the statistic $f$; this approach is also valid in the multivariate setting, as shown below. As we will see, the more substantive difficulty is that the dimension of the kernel of the $U$-statistic---namely, the subsample size $s$ of each tree---grows with the sample size.
By Proposition (ref) below, we may expand $\operatorname{RF} - \mathring \operatorname{RF}$ according to a H\"oeffding decomposition taken respect to the matrix $V^{-1}$,
where $T^{(2)}$, $T^{(3)}$, etc.\ are the second and third order projections of $T$ which obey the normal equations\footnote{ These normal equations are the main contents of Proposition (ref).}
Of course, the higher order terms $T^{(k)}$, being projections of $T$, also satisfy
These two relationships, used with (ref) and (ref), imply that
The remainder of this section centers around proving that the quantity on the RHS converge to zero. For comparison, a central result of crf (using our notation) is the bound on the diagonal elements of $\operatorname{\mathbf{Var}} \mathring T^{-1}$ and $\operatorname{\mathbf{Var}} T$. Specifically, the authors obtain
for some constant $c$. As we will see in the next section, the required bound on the trace will follow by developing bounds on the off-diagonal elements of $\operatorname{\mathbf{Var}} \mathring T$, i.e., bounds on the covariance between random forest estimates at different test points (see discussion following Proposition (ref)).
The aim of this section is to establish asymptotic bounds on the off-diagonal elements of the the covariance matrix $\operatorname{\mathbf{Var}} \mathring T$. We shall show that when the tree estimator employs splitting algorithms satisfying suitable stability conditions, we have the asymptotic behavior
Before proceeding, we first show that that this bound, coupled with control on the diagonal terms (ref), suffices to establish the trace bound in (ref) vanish.
Proposition (ref) establishes $\mathring T$ as the central object of study. Recall that $T$ is the tree estimator while $\mathring T$ is its Hajek projection; in other words, $\operatorname{\mathbf{Var}} \mathring T$ is not covariance matrix of tree estimates. However, our result will demonstrate the the asymptotic normality $V^{-1}(\operatorname{RF} - \mu)$, where $V$, the variance of the Hajek projection $\mathring\operatorname{RF}$, is given in terms of $\operatorname{\mathbf{Var}} \mathring T$ (c.f., (ref)). Therefore, (a rescaled version of) $\mathring T$ is precisely the object needed to conduct inference on the random forest. In particular, combining (ref) and (ref) yields the fact that $\operatorname{\mathbf{Var}} \mathring T$---and hence the asymptotic variance of $\operatorname{RF}$---is diagonally dominant (i.e., tending to a diagonal matrix in the limit).
We may always relabel indices so that the tree is grown on the observations $Z_1, \dots, Z_s$. To establish the bound (ref), start with the definition
To develop the term on the RHS, use the orthogonality condition for conditional expectation
Since the tree algorithm is honest, the difference $\operatorname{\mathbf{E}}(T \mid Z_1) - \operatorname{\mathbf{E}}(T \mid X_1)$ simplifies, so that for each $1 \leq k \leq q$,
where $T_k$ is the tree estimate at $x_k$, and $I_k$ the indicator for whether $X_1$ and $x_k$ belong to the same terminal node. Therefore, off-diagonal entry at $(k, l)$ of $\operatorname{\mathbf{Var}}[\operatorname{\mathbf{E}}(T \mid Z_1) - \operatorname{\mathbf{E}}(T \mid X_1)]$ is equal to
If we expand the terms in the integrand, every term will have the shape
for some multinomial $p$ of degree at most two. Since we have assumed that $\operatorname{\mathbf{E}}(Y^2 \mid X = x)$ and $\operatorname{\mathbf{E}}(Y^2 \mid X = x)$ are continuous and hence bounded, $\operatorname{\mathbf{E}}(p \mid X_1 = x)$ is also bounded. Using the Law of Iterated Expectations to evaluate (ref) then shows that it is bounded by a constant times
Recall that $I_k$ and $I_l$ are indicator variables for whether $X_1$ belong to the same hypercube as $x_k$ and $x_l$, respectively. Therefore, $\operatorname{\mathbf{E}}(I_k \mid X_1)$ is the probability that the first observation is used for the prediction at $x_k$, and likewise for $\operatorname{\mathbf{E}}(I_l \mid X_1)$. Intuitively, this only happens when $X_1$ is near $x_k$ (respectively, $x_l$): since $x_k \neq x_l$, $X_1$ cannot be near to both, meaning that the product $\operatorname{\mathbf{E}}(I_k \mid X_1)\operatorname{\mathbf{E}}(I_l \mid X_1)$ is small.
This proposition shows that the contribution of $\operatorname{\mathbf{Var}}[\operatorname{\mathbf{E}}(T \mid Z_1) - \operatorname{\mathbf{E}}(T \mid X_1)]$ to the cross covariances of $\operatorname{\mathbf{Var}} \operatorname{\mathbf{E}}(T \mid Z_1)$ is small, in particular smaller than the required bound $(\log^p s \cdot s)^{-1}$. The requirement that $\delta > 1/2$, while needed for the proof to go through, is almost certainly not needed in practice. The reason is that our proof uses $\delta > 1/2$ to derive a uniform bound on the quantity
while proposition demands a bound on its expectation. Indeed, in the extreme case $x = 0$ and $\bar x = (1, \dots, 1)^{\intercal}$, it is easy to see that the expectation satisfies the required bound even when $\delta \leq 1/2$.
Furthermore, our proof is agnostic to the exact splitting rule used by the base tree learner and uses only “random splits” (c.f., Assumption 2) in derive the required bounds. With a specific splitting rule (e.g., (ref)) and a specific data distribution, the expectation $M(x,\bar x)$ will be smaller than that predicted by (ref). In light of this, an alternative to our cyclic splitting assumption is to assume the high level condition that the splitting algorithm and data generating process confer the bound
We turn next to bound the off-diagonal terms in $\operatorname{\mathbf{Var}}[\operatorname{\mathbf{E}}(T \mid X_1)]$. As in the statement of Proposition (ref), it will be convenient to slightly change notations. We fix $1 \leq k \neq l \leq q$, and use the notation $x \mapsto x_k$, $\bar x \mapsto x_l$, with $x_1$ denoting the value of $X_1$. The goal of this section is to establish the bound
where $T$ and $\bar T$ are the tree estimates at $x$ and $\bar x$, and $\mu$ and $\bar \mu$ are the (unconditional) expectations of $T$ and $\bar T$. (Note that we have also changed the notation of $T$ slightly; before $T$ was the $q$-dimensional vector of the estimate at all points; for this section, it is the pointwise estimate at $x = x_k$.)
The quantity $\operatorname{\mathbf{E}}(T \mid X_1 = x_1) - \mu = \operatorname{\mathbf{E}}(T \mid X_1 = x_1) - \operatorname{\mathbf{E}}(T)$ measures the degree of “information” that the location of a single observation $X_1$ carries for the output of the tree at $x$. Intuitively, when $x_1$ is near $x$, the effect of $x_1$ on the leaf node containing $x$ is more pronounced and we expect $\operatorname{\mathbf{E}}(T \mid X_1 = x_1) - \mu \approx \operatorname{\mathbf{E}}(Y \mid X_1 = x_1) - \mu$. Conversely, when $x_1$ is far from $x$, then its effect in determining the location of the leaf node containing $x$ diminishes, and $\operatorname{\mathbf{E}}(T \mid X_1 = x_1) - \mu \approx 0$.
The key in making the above intuition precise is to keep track of the $X_1 = x_1$ leaves the intermediate partition containing $x$, where “intermediate partition” refers to the nodes created during the splitting process that are not necessarily terminal. We will see that once $X_1$ and $x$ are separated, its effect on the prediction decreases.
Towards this end, fix $x$ and let $\Pi$ denote the terminal node containing $x$; $\Pi$ is a random subset of $\mathcal{X}$ created by axis aligned splits. By Assumption 4, the set of potential splits does not depend on the sample (in particular, it does not depend on $X_1$). Moreover, splitting ceases after no more than $s$ times, regardless of the subsample $X_1, \dots, X_s$, as each split reduces the number of observations in its two children nodes by at least one. Therefore, $\Pi$ takes on only finitely many possible values, and we may write
where $\mu_{\pi} = \operatorname{\mathbf{E}}(T \mid \Pi = \pi)$ and $\mu_{\pi}' = \operatorname{\mathbf{E}}(T \mid \Pi = \pi, X_1 = x)$.
The hyperrectangle $\Pi$ is determined by the recursive splitting procedure used to grow the tree, and there is a natural correspondence between (ref) and a certain “expectation” taken over the directed acyclic graph (DAG) defined in the following way. Let $[0, 1]^p$ be the root of the DAG; for every potential split at $[0,1]^p$, there is a directed edge to a new vertex, where the vertex is the hyperrectangle that contains $x$. If the node represented by a vertex is one of the possible values of $\Pi$, then that vertex is a leaf in the DAG and has no outgoing edges; other vertices carry an outgoing edge for each potential split at that node, with each edge going to another vertex which is again a hyperrectangle containing $x$.
The previous definition determines the DAG recursively: each vertex in the DAG is a node containing $x$, with terminal vertices corresponding to terminal nodes. To each terminal vertex $v$, we associate the value $f(v) \coloneq \mu_{\pi}$ as in (ref). In addition, each edge $e = (v \to w)$ corresponds to a split at a node $v$ producing a halfspace $w$ of $v$; associate with this edge a “transition probability”
Given the transition probabilities, the value $f$ may be extended to each vertex $v$ recursively via the formula
We refer to $f$ as the continuation value at $v$, and by construction we have
Alternatively, if we had assigned the values $f'(v) = \mu'_v$ to each terminal vertex and used the transition probabilities
then we recover $\operatorname{\mathbf{E}}(T \mid X_1 = x_1) = f'([0,1]^p)$ after extending $f'$ in the same way as $f$. In other words, bounding $\operatorname{\mathbf{E}}(T \mid X_1 = x_1) - \operatorname{\mathbf{E}}(T)$ requires bounding the difference between the two types of continuation values.
We will need to assume that $p'(e) \approx p(e)$; that is, conditioning on a single observation will not affect the probability that a particular split is chosen (i.e., that a particular split is optimal at its node). This is a natural assumption in that the optimal split is computed using all the observations in a particular node, so that conditioning on a single observation should have relatively little effect. This will depend on the specifics on the splitting algorithm used to a construct a tree, and rules which satisfying the following assumption to be “stable splitting rules.” \begin{splitting-stability} For any node $v$, the total variation distance between the distributions $\{ p(e) \}_{e: v \to w}$ and $\{ p'(e) \}_{e: v \to w}$ is bounded by the volume of $v$. Specifically, there exists some $\epsilon>0$ such that for all $v$,
Here, $|v|$ denotes the volume of the hyperrectangle at $v$, i.e.,
\end{splitting-stability}
Since the distribution of $X$ has a density that is bounded above and below, a simple H\"oeffding bound shows that the number of sample points in $v$ is bounded above and below by $s|v|$, with the constants adjusted so that the failure probability is less than\footnote{For example, the probability that a binomial random variable $B(n,p)$ deviates from $np$ by more than $\sqrt{n \log n}$ is less than $C/s^2$ for some constant $C$.} $1/s^2$. Since this is smaller than the required bound $(\log^p s \cdot s)^{-1}$, we may interpret $s|v|$ in (ref) to be the number of samples in $|v|$ without loss of generality. Relatedly, Wager2015AdaptiveCO's Lemma 12 (see also the proof of Lemma 2 in crf) extends to this fact to be uniform across nodes.
The stability assumption places a restriction on procedure used to select optimal splits: namely, if the decision is made on the basis of $m$ points, then conditioning on any one of the points changes the optimal split with probability bounded by $m^{-(1+\epsilon)}$. In practice, most splitting procedures satisfy a stronger bound. A set of sufficient conditions is given in the following proposition.
In general, the conditions in Proposition (ref) are sufficient to guarantee an exponential bound instead of a polynomial one as in (ref). Thus, Proposition (ref) should be viewed as simply as providing a plausibility argument that stable splitting rules are commonly encountered in practice.
The next proposition shows that bounds on splitting probabilities automatically imply a related bound on the continuation values.
The splitting stability assumption stipulates that $\Delta(z) = z^{-(1+\epsilon)}$, where the factor $\epsilon$ allows us to ignore the extra logarithm. In that case, we may put the bounds on $\operatorname{TV}(p, p')$ and $|f-f'|$ together and establish required bound on $\operatorname{\mathbf{Var}} \operatorname{\mathbf{E}}(T \mid X_1)$.
Recall that Proposition (ref) required that $\delta > 1/2$. Since $\alpha < 1/2$ by definition, the requirement that $\delta > 1-\alpha$ in Proposition (ref) is more restrictive. Just like Proposition (ref), we argue that this requirement is plausibly looser in applications. The reason is that it is used to give the following bound on hyperrectangles $v$ created after $L$ splits
The RHS appears since potential split may reduce the volume of a node by at most $\alpha$: but only an exponentially small (i.e., $2^{-L}$) proportion of nodes are the result of taking the smallest possible split $L$ times! The “average” node has volume approximately $(1/2)^L$, so that $\delta > 1/2$ may be more appropriate.
Combining Propositions (ref) and (ref) with equations (ref) and (ref) yields the desired bound (ref) on the off-diagonal terms of $\operatorname{\mathbf{Var}} \mathring T$ discussed at the beginning of this section. Therefore, Proposition (ref) applies, and the joint normality of the random forest estimator is established.
The previous sections focused on deriving the asymptotic normality result
Recall our standing convention that $\operatorname{RF}\in \mathbf{R}^q$ is the random forest estimate at $x_1, \dots, x_q$ and $\mu$ is its expectation. According to (ref), the target function of the random forest is actually is $m(x) = \operatorname{\mathbf{E}}(Y \mid X = x)$. The results in crf show that $(\operatorname{RF}(x) - m(x))/\sqrt{V} \xRightarrow{\;\rm dist\;} N(0, 1)$ pointwise for each $x \in \mathcal{X}$; since we have shown that $V$ is diagonally dominant in that its off-diagonal terms vanish relative to the diagonal, the pointwise result carries over to our multivariate setting, and $V^{-1/2}(\operatorname{RF} - m) \xRightarrow{\;\rm dist\;} N(0, I)$, where $m = (m(x_1), \dots, m(x_q))$.
Moreover, crf proposes a jackknife estimator that can consistently estimate $\sqrt{V}$ in the scalar case. Our diagonal dominance result implies that the random forest estimates at $x$ and $\bar x \in [0, 1]^p$ are independent in the limit $n\to\infty$,
so that the jackknife estimator for the scalar case may be fruitfully applied to obtain confidence bands for functionals of the random forest estimates (i.e., expressions involving the estimate at more than one point).
The accuracy of the approximation above above depends on the decay of the off-diagonal terms in finite samples. In this section, we provide a “back of the envelope” bound for the covariance term that may be useful for practitioners. We stress that the following calculations are (mostly) heuristics: as we have shown above, the covariance term depends on quantities such as $M(x, \bar x)$, which is in turn heavily dependent on the exact mechanics of the underlying splitting algorithm. Since our aim is to produce a “usable” result, we will dispense with rigorous analysis in the remainder of this section.
To begin, the proofs of Propositions (ref) and (ref) showed that the asymptotic variance $V$ has off-diagonal terms which are upper bounded by\footnote{This is a very crude upper bound as we have dropped the quantity $\Delta(\alpha^{\ell} s)$ from the infinite series.}
where $p_{\ell} = \operatorname{\mathbf{P}}(L \geq \ell)$ is the probability that $x$ and $x_1$ are not separated after $\ell$ splits and likewise for $\bar p_{\ell} = \operatorname{\mathbf{P}}(\bar L \geq \ell)$. That is, $L$ is the number of splits before which $x$ and $x_1$ belong to the same partition. If we denote by $I$ (resp.\ $\bar I$) the indicator variable that $X_1$ is in the terminal node of $x$ (resp.\ $\bar x$), then the events $\{ I=1 \}$ and $\{L = \log_2 s\}$ are equal, so that
Replacing the inequality with an approximation, we have $\operatorname{\mathbf{E}} L = (\log s) \operatorname{\mathbf{E}}(I \mid X_1 = x_1)$. All of this shows that the covariance term is bounded by
Towards a useful heuristic, we will consider a bound on the correlation instead of the covariance. In our notation, the result of crf lower bounds $M(x, x)$ (and $M(\bar x, \bar x)$), while our paper provides an upper bound on $M(x, \bar x)$. Ignoring the logarithmic terms, we have
Recall that $M(x, \bar x) = \operatorname{\mathbf{E}}[\operatorname{\mathbf{E}}(I \mid X_1) \operatorname{\mathbf{E}}(\bar I \mid X_1)]$, which decays as $\bar x$ moves away from $x$. Using the previous expression (note that $M(x,x) \approx M(\bar x, \bar x)$ due to symmetry between $x$ and $\bar x$), we can bound the correlation from purely geometric considerations. Since the integrand
decays as $X_1$ moves away from $x$ (and $\bar x$), we may imagine that its integral
has the largest contributions for points $x_1$ near $x$, say those points in a $L_{\infty}$-box of side lengths $d$ with volume $d^p$, i.e., $\{ y \in [0,1]^p : \| x - y \| \leq d/2\}$. If we accept this, then the contributions for the integral $M(x, \bar x)$ would come from points that are with $d/2$ of both $x$ and $\bar x$, and to a first degree approximation, the volume of these points $\{ y \in [0,1]^p : \| x - y \|_{\infty} \leq d/2, \| \bar x - y \| \leq d/2 \}$ is
where the approximation is accurate if $|z_i| \ll 1$. Dividing through by $d^p$, the proportion of the volume of the latter set is $1 - \frac{1}{d}\| x - \bar x \|_1$, which leads to the heuristic
The RHS has the correct scaling when $x = \bar x$, i.e., the correlation equals one when $\| x - \bar x \|_1 = 0$. To maintain correct scaling at the other extreme with $\| x - \bar x \|_1 = p$, we should take $c = 1/p$, so that
Of course, this heuristic must be incorrect in that it does not depend on $s$; our theoretical results show that even for non-diametrically opposed points, the correlation drops to zero as $s \to \infty$. Therefore, another recommendation is to use
for some $\epsilon > 0$, where dependence $s^{\epsilon}$ comes from considering the decay of $M(x, \bar x)$ as $\bar x$ moves away from $x$ (c.f.\ the proof Proposition (ref)).
In this section, we discuss the results of simulations calculating the correlation structure. For our experiment, we set $p = 2$, so that the covariates $X$ are distributed on the unit square. The distribution of $X$ is chosen to be “four-modal”
and $\bar N$ denotes a truncated multivariate Gaussian distribution on the unit square.\footnote{ That is, $\bar N(\mu, \Sigma)$ denotes the conditional distribution of $x \sim N(\mu, \Sigma)$ on the event $x \in [0, 1]^2$.} Thus, $X$ has a bounded density on the unit square, and has four peaks at $\mu_1, \dots, \mu_4$. The distribution of $Y$ conditional on $X = (x_1, x_2)$ is
The random splitting probability is $\delta = 1/2$, and the regularity parameters are $\alpha = 0.01$ and $k = 1$, so that the tree is grown to the fullest extent (i.e., terminal nodes may contain a single observation), with each terminal node lying on the $101 \times 101$ grid of the unit square. For each sample size $n$, five thousand trees are grown, and the estimates are aggregated to compute the correlation.
Figure (ref) plots the correlation of estimates at $x$ and $\bar x$ as a function of the $L_1$ norm $\|x - \bar x\|_1$. The calculation is performed by first fixing $x$, then calculating the sample correlation (across five thousand trees) as $\bar x$ ranges over each cell: the correlation is associated with the $L_1$ norm $\|x - \bar x\|_1$. This process is then repeated by varying the reference point $x$, and the correlation at $\| x - \bar x \|_1$ is the average of the correlations observed.
The figure demonstrates that the linear heuristic (ref) given in the previous section is conservative: it is evident that correlation decreases super-linearly as $x$ and $\bar x$ become separated.
Figure (ref) plots the correlation on a logarithmic scale, which shows that that correlation decay is exponential in a neighborhood of unity. In other words, simulations suggest that the correct heuristic is of the shape
for a suitable $\lambda$.
Random forests and tree-based methods form an important part of of an applied data analysis toolkit. In this paper, we studying the covariance between random forest estimates at several points. We develop a novel construction of a directed acyclic graph that keeps track of the splitting probabilities when knowledge of one point is known (Propositions (ref) and (ref)). As part of the proof, we establish stability properties of a class of splitting rules (see Proposition (ref)). We also identify (Proposition (ref)) $M(x, \bar x)$, which (roughly) captures the likelihood of two points belonging to the same terminal node, as a key quantity in controlling the off-diagonal term of the covariance matrix of the multivariate random forest.
In this way, this paper provides the a theoretical basis for performing inferences on functionals of target function (e.g., a heterogeneous treatment effect) when the functional is based on values of the target function at multiple points in the feature space. Specifically, we show that the covariance vanishes in the limit relative to the variance, and provide heuristics on the size of the correlation in finite samples.
We close with discussing a couple avenues for future research. The first is extending our framework to cover categorical or discrete-valued features. Here, new assumptions would be required in order to maintain the guarantee that node sizes are “not too small.” Second, our bounds---after potential improvements---on the covariance matrix of the random forest may be used with the recent results of chernozhukov2017,chen2018 in order to provide finite sample Gaussian approximations. This would provide a sounder theoretical underpinning of for our heuristics, and increase the usefulness of this paper for practitioners.