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.
151,341 characters · 21 sections · 157 citation commands
Generalized Random Forests
Random forests, introduced by breiman2001random, are a widely used algorithm for statistical learning. Statisticians usually study random forests as a practical method for non-parametric conditional mean estimation: Given a data-generating distribution for $(X_i, \, Y_i) \in \xx \times \RR$, forests are used to estimate $\mu(x) = \EE{Y_i \cond X_i = x}$. Several theoretical results are available on the asymptotic behavior of such forest-based estimates \smash{$\hmu(x)$}, including consistency arlot2014analysis,biau2008consistency,biau2012analysis,denil2014narrowing,lin2006random, scornet2015consistency,wager2015uniform, second-order asymptotics mentch2016quantifying, and confidence intervals wager2015estimation.
This paper extends Breiman's random forests into a flexible method for estimating any quantity $\theta(x)$ identified via local moment conditions. Specifically, given data $(X_i, \, O_i) \in \xx \times \oo$, we seek forest-based estimates of $\theta(x)$ defined by a local estimating equation of the form
\sloppy{where $\psi(\cdot)$ is some scoring function and $\nu(x)$ is an optional nuisance parameter. This setup encompasses several key statistical problems. For example, if we model the distribution of $O_i$ conditionally on $X_i$ as having a density $f_{\theta(x), \, \nu(x)}(\cdot)$ then, under standard regularity conditions, the moment condition (ref) with \smash{$\psi_{\theta(x), \, \nu(x)}(O) = \nabla \log\p{f_{\theta(x), \, \nu(x)}(O_i)}$} identifies the local maximum likelihood parameters $(\theta(x), \, \nu(x))$. More generally, we can use moment conditions of the form (ref) to identify conditional means, quantiles, average partial effects, etc., and to develop robust regression procedures via Huberization. Our main substantive application of generalized random forests involves heterogeneous treatment effect estimation with instrumental variables.}
Our aim is to build a family of non-parametric estimators that inherit the desirable empirical properties of regression forests---such as stability, ease of use, and flexible adaptation to different functional forms as in, e.g., biau2016random or varian2014big---but can be used in the wide range of statistical settings characterized by (ref) in addition to standard conditional mean estimation. This paper addresses the resulting conceptual and methodological challenges and establishes formal asymptotic results.
Regression forests are typically understood as ensemble methods, i.e., forest predictions \smash{$\hmu(x)$} are written as the average of $B$ noisy tree-based predictors \smash{$\hmu_b(x)$}, \smash{$\hmu(x) = B^{-1} \sum_{b = 1}^B \hmu_b(x)$}; and, because individual trees \smash{$\hmu_b(x)$} have low bias but high variance, such averaging meaningfully stabilizes predictions buhlmann2002analyzing,scornet2015consistency. However, noisy solutions to moment equations as in (ref) are generally biased, and averaging would do nothing to alleviate the bias.
To avoid this issue, we cast forests as a type of adaptive locally weighted estimators that first use a forest to calculate a weighted set of neighbors for each test point $x$, and then solve a plug-in version of the estimating equation (ref) using these neighbors. Section (ref) gives a detailed treatment of this perspective. This locally weighting view of random forests was previously advocated by hothorn2004bagging in the context of survival analysis and by meinshausen2006quantile for quantile regression, and also underlies theoretical analyses of regression forests lin2006random. For conditional mean estimation, the averaging and weighting views of forests are equivalent; however, once we move to more general settings, the weighting-based perspective proves substantially more effective, and also brings forests closer to the literature on local maximum likelihood estimation fan1996local,loader1999local,newey1994kernel,stone1977consistent,tibshirani1987local.
A second challenge in generalizing forest-based methods is that their success hinges on whether the adaptive neighborhood function obtained via partitioning adequately captures the heterogeneity in the underlying function $\theta(x)$ we want to estimate. Even within the same class of statistical tasks, different types of questions can require different neighborhood functions. For example, suppose that two scientists are studying the effects of a new medical treatment: One is looking at how the treatment affects long-term survival, and the other at its effect on the length of hospital stays. It is plausible that the treatment heterogeneity in each setting would be based on disparate covariates, e.g., a patient's smoking habits for long-term survival, and the location and size of the hospital for the length of stay.
Thus, each time we apply random forests to a new scientific task, it is important to use rules for recursive partitioning that are able to detect and highlight heterogeneity in the signal the researcher is interested in. In prior work, such problem-specific rules have largely been designed on a case by case basis. Although the CART rules of breiman1984classification have long been popular for classification and regression tasks, there has been a steady stream of papers proposing new splitting rules for other problems, including athey2016recursive and su2009subgroup for treatment effect estimation, beygelzimer2009offset and kallus2016learning for personalized policy allocation, and gordon1985tree, leblanc1992relative, molinaro2004tree as well as several others for survival analysis. zeileis2008model propose a method for constructing a single tree for general maximum likelihood problems, where splitting is based on hypothesis tests for improvements in goodness of fit.
In contrast, we seek a unified, general framework for computationally efficient problem-specific splitting rules, optimized for the primary objective of capturing heterogeneity in a key parameter of interest. In the spirit of gradient boosting friedman2001greedy, our recursive partitioning method begins by computing a linear, gradient-based approximation to the non-linear estimating equation we are trying to solve, and then uses this approximation to specify the tree-split point. Algorithmically, our procedure reduces to iteratively applying a labeling step where we generate pseudo-outcomes by computing gradients using parameters estimated in the parent node, and a regression step where we pass this labeled data to a standard CART regression routine. Thus, we can make use of pre-existing, optimized tree software to execute the regression step, and obtain high quality neighborhood functions while only using computational resources comparable to those required by standard CART algorithms. In line with this approach, our generalized random forest software package builds on the carefully optimized ranger implementation of regression forest splitting rules wright2015ranger.
Moment conditions of the form (ref) typically arise in scientific applications where rigorous statistical inference is required. The bulk of this paper is devoted to a theoretical analysis of generalized random forests, and to establishing asymptotic consistency and Gaussianity of the resulting estimates \smash{$\htheta(x)$}. We also develop methodology for asymptotic confidence intervals. Our analysis is motivated by classical results for local estimating equations, in particular newey1994kernel, paired with machinery from wager2015estimation to address the adaptivity of the random forest weighting function.
The resulting framework presents a flexible method for non-parametric statistical estimation and inference with formal asymptotic guarantees. In this paper, we develop applications to quantile regression, conditional average partial effect estimation and heterogeneous treatment effect estimation with instrumental variables; however, there are many other popular statistical models that fit directly into our framework, including panel regression, Huberized robust regression, models of consumer choice, etc. In order to fit any of these models with generalized random forests, the analyst simply needs to provide the problem-specific routines to calculate gradients of the moment conditions evaluated at different observations in the dataset for the “label” step of our algorithm. Moreover, we emphasize that our method is in fact a proper generalization of regression forests: If we apply our framework to build a forest-based method for local least-squares regression, we exactly recover a regression forest. A high-performance software implementation of generalized random forests, grf for R and C++, is available from CRAN.
The idea of local maximum likelihood (and local generalized method of moments) estimation has a long history, including fan1998local, newey1994kernel, staniswalis1989kernel, stone1977consistent, tibshirani1987local and lewbel2007local. In economics, popular applications of these techniques include multinomial choice modeling in a panel data setting honore2000panel and instrumental variables regression su2013local. The core idea is that when estimating parameters at a particular value of covariates, a kernel weighting function is used to place more weight on nearby observations in the covariate space. A challenge facing this approach is that if the covariate space has more than two or three dimensions, performance can suffer due to the “curse of dimensionality” robins1997toward.
Our paper replaces the kernel weighting with forest-based weights, that is, weights derived from the fraction of trees in which an observation appears in the same leaf as the target value of the covariate vector. The original random forest algorithm for non-parametric classification and regression was proposed by breiman2001random, building on insights from the ensemble learning literature amit1997shape,breiman1996bagging,dietterich2000experimental,ho1998random. The perspective we take on random forests as a form of adaptive nearest neighbor estimation, however, most closely builds on the proposals of hothorn2004bagging and meinshausen2006quantile for forest-based survival analysis and quantile regression. This adaptive nearest neighbors perspective also underlies several statistical analyses of random forests, including arlot2014analysis, biau2010layered, and lin2006random.
Our gradient-based splitting scheme draws heavily from a long tradition in the statistics and econometrics literatures of using gradient-based test statistics to detect change points in likelihood models andrews1993tests,hansen1992testing,hjort2002tests,nyblom1989testing, ploberger1992cusum,zeileis2005unified,zeileis2007generalized. In particular, zeileis2008model consider the use of such methods for model-based recursive partitioning. Our problem setting differs from the above in that we are not focused on running a hypothesis test, but rather seek an adaptive nearest neighbor weighting that is as sensitive as possible to heterogeneity in our parameter of interest; we then rely on the random forest resampling mechanism to achieve statistical stability mentch2016quantifying,scornet2015consistency,wager2015estimation. In this sense, our approach is related to gradient boosting friedman2001greedy, which uses gradient-based approximations to guide a greedy, non-parametric regression procedure.
Our asymptotic theory relates to an extensive recent literature on the statistics of random forests, most of which focuses on the regression case arlot2014analysis,biau2012analysis,biau2008consistency,biau2016random, buhlmann2002analyzing,denil2014narrowing,geurts2006extremely, ishwaran2010consistency,lin2006random,meinshausen2006quantile,mentch2016quantifying, scornet2015consistency,sexton2009standard,wager2015estimation,wager2015uniform,zhu2015reinforcement. Our present paper complements this body of work by showing how methods developed to study regression forests can also be used understand estimated solutions to local moment equations obtained via generalized random forests.
In standard classification or regression forests as proposed by breiman2001random, the prediction for a particular test point $x$ is determined by averaging predictions across an ensemble of different trees amit1997shape,breiman1996bagging,dietterich2000experimental,ho1998random. Individual trees are grown by greedy recursive partitioning, i.e., we recursively add axis-aligned splits to the tree, where each split it chosen to maximize the improvement to model fit breiman1984classification; see Figure (ref) in the Appendix for an example of a tree. The trees are randomized using bootstrap (or subsample) aggregation, whereby each tree is grown on a different random subset of the training data, and random split selection that restricts the variables available at each step of the algorithm. For an introductory overview of random forests, we recommend the chapter of hastie2009elements dedicated to the method. As discussed below, in generalizing random forests, we preserve several core elements of Breiman's forests---including recursive partitioning, subsampling, and random split selection---but we abandon the idea that our final estimate is obtained by averaging estimates from each member of an ensemble. Treating forests as a type of adaptive nearest neighbor estimator is much more amenable to statistical extensions.
Suppose that we have $n$ independent and identically distributed samples, indexed $ i = 1, \, ..., \, n$. For each sample, we have access to an observable quantity $O_i$ that encodes information relevant to estimating $\theta(\cdot)$, along with a set of auxiliary covariates $X_i$. In the case of non-parametric regression, this observable just consists of an outcome $O_i = \cb{Y_i}$ with $Y_i \in \RR$; in general, it may contain richer information. In the case of treatment effect estimation with exogenous treatment assignment, $O_i = \cb{Y_i, \, W_i}$ also includes the treatment assignment $W_i$. Given this type of data, our goal is to estimate solutions to local estimation equations of the form \smash{$\mathbb{E}[\psi_{\theta(x), \, \nu(x)}\p{O_i} \cond X_i = x] = 0$} for all $\in \xx$, where $\theta(x)$ is the parameter we care about and $\nu(x)$ is an optional nuisance parameter.
One approach to estimating such functions $\theta(x)$ is to first define some kind of similarity weights $\alpha_i(x)$ that measure the relevance of the $i$-th training example to fitting $\theta(\cdot)$ at $x$, and then fit the target of interest via an empirical version of the estimating equation fan1998local,newey1994kernel,staniswalis1989kernel,stone1977consistent,tibshirani1987local:
When the above expression has a unique root, we can simply say that \smash{$(\htheta(x), \, \hnu(x))$} solves \smash{$\sum_{i = 1}^n \alpha_i(x) \, \psi_{\htheta(x), \, \hnu(x)}\p{O_i} = 0$}. The weights $\alpha_i(x)$ used to specify the above solution to the heterogeneous estimating equation are traditionally obtained via a deterministic kernel function, perhaps with an adaptively chosen bandwidth parameter hastie2009elements. Although methods of the above kind often work well in low dimensions, they are sensitive to the curse of dimensionality.
Here, we seek to use forest-based algorithms to adaptively learn better, problem-specific, weights $\alpha_i(x)$ that can be used in conjunction with (ref). As in hothorn2004bagging and meinshausen2006quantile, our generalized random forests obtain such weights by averaging neighborhoods implicitly produced by different trees. First, we grow a set of $B$ trees indexed by $b = 1, \, ..., \, B$ and, for each such tree, define $L_b(x)$ as the set of training examples falling in the same “leaf” as $x$. The weights $\alpha_i(x)$ then capture the frequency with which the $i$-th training example falls into the same leaf as $x$:
These weights sum to 1, and define the forest-based adaptive neighborhood of $x$; see Figure (ref) for an illustration of this weighting function.
There are some subtleties in how the sets $L_b(x)$ are defined---in particular, as discussed in Section (ref), our construction will rely on both subsampling and a specific form of sample-splitting to achieve consistency---but at a high level the estimates \smash{$\htheta(x)$} produced by a generalized random forests are simply obtained by solving (ref) with weights (ref).
Finally, for the special case of regression trees, our weighting-based definition of a random forest is equivalent to the standard “average of trees” perspective taken in breiman2001random: If we estimate the conditional mean function $\mu(x) = \EE{Y_i \cond X_i = x}$, as identified in (ref) using $\psi_{\mu(x)}(Y_i) = Y_i - \mu(x)$, then we see that \smash{$\sum_{i = 1}^n \frac{1}{B} \sum_{b = 1}^B \alpha_{bi}(x) \p{Y_i - \hmu(x)} = 0$} if and only if \smash{$\hmu(x) = \frac{1}{B} \sum_{b = 1}^B \hmu_b(x)$}, where \smash{$\hmu_b(x) = \sum_{\cb{i : X_i \in L_b(x)}} Y_i \, \big/ \, \abs{L_b(x)}$} is the prediction made by a single CART regression tree.
We seek trees that, when combined into a forest, induce weights $\alpha_i(x)$ that lead to good estimates of $\theta(x)$. The main difference between random forests relative to other non-parametric regression techniques is their use of recursive partitioning on subsamples to generate these weights $\alpha_i(x)$. Motivated by the empirical success of regression forests across several application areas, our approach mimics the algorithm of breiman2001random as closely as possible, while tailoring our splitting scheme to focus on heterogeneity in the target functional $\theta(x)$.
Just like in Breiman's forests, our search for good splits proceeds greedily, i.e., we seek splits that immediately improve the quality of the tree fit as much as possible. Every split starts with a parent node $P \subseteq \xx$; given a sample of data $\mathcal{J}$, we define \smash{$(\htheta_P, \, \hnu_P)(\mathcal{J})$} as the solution to the estimating equation, as follows (we suppress dependence on $\mathcal{J}$ when unambiguous):
We would like to divide $P$ into two children $C_1, \, C_2 \subseteq \xx$ using an axis-aligned cut such as to improve the accuracy of our $\theta$-estimates as much as possible; formally, we seek to minimize \smash{$\err\p{C_1, \, C_2}$} defined as \smash{$\err\p{C_1, \, C_2} = \sum_{j = 1, \, 2} \mathbb{P}[X \in C_j \cond X \in P] \, \mathbb{E}[(\htheta_{C_j}(\mathcal{J}) - \theta(X))^2 \cond X \in C_j]$}, where \smash{$\htheta_{C_j}(\mathcal{J})$} are fit over children $C_j$ in analogy to (ref), and expectations are taken over both the randomness in \smash{$\htheta_{C_j}(\mathcal{J})$} and a new test point $X$.
Many standard regression tree implementations, such as CART breiman1984classification, choose their splits by simply minimizing the in-sample prediction error of the node, which corresponds to \smash{$\err\p{C_1, \, C_2}$} with plug-in estimators from the training sample. In the case of estimating the effect of a binary treatment, athey2016recursive study sample-splitting trees, and propose an unbiased, model-free estimate of \smash{$\err\p{C_1, \, C_2}$} using an overfitting penalty in the spirit of mallows1973some. In our setting, however, this kind of direct loss minimization is not an option: If $\theta(x)$ is only identified through a moment condition, then we do not in general have access to unbiased, model-free estimates of the criterion \smash{$\err\p{C_1, \, C_2}$}. To address this issue, we rely on the following more abstract characterization of our target criterion.
Motivated by this observation, we consider splits that make the above $\Delta$-criterion (ref) large. A special case of the above idea also underlies the splitting rule for treatment effect estimation proposed by athey2016recursive. At a high level, we can think of this $\Delta$-criterion as favoring splits that increase the heterogeneity of the in-sample $\theta$-estimates as fast as possible. The dominant bias term in \smash{$\err\p{C_1, \, C_2}$} is due to the sampling variance of regression trees, and is the same term that appears in the analysis of athey2016recursive. Including this error term in the splitting criterion may stabilize the construction of the tree, and further it can prevent the splitting criterion from favoring splits that make the model difficult to estimate.
The above discussion provides conceptual guidance on how to pick good splits. But actually optimizing the criterion $\Delta(C_1, \, C_2)$ over all possible axis-aligned splits while explicitly solving for $\htheta_{C_1}$ and $\htheta_{C_2}$ in each candidate child using an analogue to (ref) may be quite expensive computationally. To avoid this issue, we instead optimize an approximate criterion \smash{$\tDelta(C_1, \, C_2)$} built using gradient-based approximations for $\htheta_{C_1}$ and $\htheta_{C_2}$. For each child $C$, we use $\ttheta_C \approx \htheta_C$ as follows: We first compute $A_P$ as any consistent estimate for the gradient of the expectation of the $\psi$-function, i.e., $\nabla \mathbb{E}[\psi_{\htheta_P, \, \hnu_P}(O_i) \cond X_i \in P]$, and then set
where $\htheta_P$ and $\hnu_P$ are obtained by solving (ref) once in the parent node, and $\xi$ is a vector that picks out the $\theta$-coordinate from the $(\theta, \, \nu)$ vector. When the $\psi$-function itself is continuously differentiable, we use
and the quantity \smash{$\xi^\top A^{-1}_P \ \psi_{\htheta_P, \, \hnu_p}\p{O_i}$} corresponds to the influence function of the $i$-th observation for computing \smash{$\htheta_P$} in the parent. Cases where $\psi$ is non-differentiable, e.g., with quantile regression, require more care.
Algorithmically, our recursive partitioning scheme now reduces to alternatively applying the following two steps. First, in a {\bf labeling step}, we compute $\htheta_P$, $\hnu_P$, and the derivative matrix $A_P^{-1}$ on the parent data as in (ref), and use them to get pseudo-outcomes
Next, in a {\bf regression step}, we run a standard CART regression split on the pseudo-outcomes $\rho_i$. Specifically, we split $P$ into two axis-aligned children $C_1$ and $C_2$ such as to maximize the criterion
Once we have executed the regression step, we relabel observations in each child by solving the estimating equation, and continue on recursively.
For intuition, it is helpful to examine the simplest case of least-squares regression, i.e., with $\psi_{\theta(x)}(Y) = Y - \theta(x)$. Here, the labeling step (ref) doesn't change anything---we get $\rho_i = Y_i - \bY_p$, where $\bY_p$ is the mean outcome in the parent---while the second step maximizing (ref) corresponds to the usual way of making splits as in breiman2001random. Thus, the special structure of the type of problem we are trying to solve is encoded in (ref), while the second scanning step is a universal step shared across all different types of forests.
We expect this approach to provide more consistent computational performance than optimizing (ref) at each split directly. When growing a tree, the computation is typically dominated by the split-selection step, and so it is critical for this step to be implemented as efficiently as possible (conversely, the labeling step (ref) is only solved once per node, and so is less performance sensitive). From this perspective, using a regression splitting criterion as in (ref) is very desirable, as it is possible to evaluate all possible split points along a given feature with only a single pass over the data in the parent node (by representing the criterion in terms of cumulative sums). In contrast, directly optimizing the original criterion (ref) may require solving intricate optimization problems for each possible candidate split.
This type of gradient-based approximation also underlies other popular statistical algorithms, including gradient boosting friedman2001greedy and the model-based recursive partitioning algorithm of zeileis2008model. Conceptually, tree splitting bears some connection to change-point detection if we imagine tree splits as occurring at detected change-points in $\theta(x)$; and, from this perspective, our approach is closely related to standard techniques for moment-based change-point detection andrews1993tests,hansen1992testing,hjort2002tests,nyblom1989testing, ploberger1992cusum,zeileis2005unified,zeileis2007generalized.
In our context, we can verify that the error from using the approximate criterion (ref) instead of the exact $\Delta$-criterion (ref) is within the tolerance used to motivate the $\Delta$-criterion in Proposition (ref), thus suggesting that our use of (ref) to guide splitting may not result in too much inefficiency. Note that consistent estimates of $A_P$ can in general be derived directly via, e.g., (ref), without relying on Proposition (ref).
Now, given a practical splitting scheme for growing individual trees, we want to grow a forest that allows for consistent estimation of $\theta(x)$ using (ref) paired with the forest weights (ref). We expect each tree to provide small, relevant neighborhoods for $x$ that give us noisy estimates of $\theta(x)$; then, we may hope that forest-based aggregation will provide a single larger but still relevant neighborhood for $x$ that yields stable estimates $\htheta(x)$.
To ensure good statistical behavior, we rely on two conceptual ideas that have proven to be successful in the literature on forest-based least-squares regression: Training trees on subsamples of the training data mentch2016quantifying,scornet2015consistency,wager2015estimation, and a sub-sample splitting technique that we call honesty biau2012analysis,denil2014narrowing,wager2015estimation. Our final algorithm for forest-based solutions to heterogeneous estimating equations is given as Algorithm (ref); we refer to Section 2.4 of wager2015estimation for a more in-depth discussion of honesty in the context of forests. As shown in Section (ref), assuming regularity conditions, the estimates $\htheta(x)$ obtained using a generalized random forest as described in Algorithm (ref) are consistent for $\theta(x)$. Moreover, given appropriate subsampling rates, we establish asymptotic normality of the resulting forest estimates $\htheta(x)$.
We now turn to a formal characterization of generalized random forests, with the aim of establishing asymptotic Gaussianity of the \smash{$\htheta(x)$}, and of providing tools for statistical inference about $\theta(x)$. We first list assumptions underlying our theoretical results. Throughout, the covariate space and the parameter space are both subsets of Euclidean space; specifically, $\xx = [0, \, 1]^p$ and $(\theta, \, \nu) \in \bb \subset \RR^k$ for some $p, \, k > 0$, where $\bb$ is a compact subset of $\RR^k$. Moreover, we assume that the features $X$ have a density that is bounded away from 0 and $\infty$; as argued in, e.g., wager2015uniform, this is equivalent to imposing a weak dependence condition on the individual features $(X_i)_j$ because trees and forests are invariant to monotone rescaling of the features. All proofs are in the appendix.
Some practically interesting cases, such as quantile regression, involve discontinuous score functions $\psi$, which makes the analysis more intricate. Here, we follow standard practice, and assume that the expected score function,
varies smoothly in the parameters, even though $\psi$ itself may be discontinuous. For example, with quantile regression \smash{$\psi_\theta(Y) = 1\p{\cb{Y > \theta}} - (1 - q)$} is discontinuous in $q$, but \smash{$M_\theta(x) = \PP{Y > \theta \cond X = x} - (1 - q)$} will be smooth whenever $Y\cond X = x$ has a smooth density.
Our next two assumptions control regularity properties of the $\psi$-function itself. Assumption (ref) holds trivially when $\psi$ itself is Lipschitz in $(\theta, \, \nu)$ (in fact, having $\psi$ be 0.5-H{\"o}lder would be enough), while Assumption (ref) is used to show that a certain empirical process is Donsker. Examples are given at the end of this section.
All the previous assumptions only deal with local properties of the estimating equation, and can be used to control the behavior of \smash{$(\htheta(x), \, \hnu(x))$} in a small neighborhood of the population parameter value $(\theta(x), \, \nu(x))$. Now, to make any use of these assumptions, we first need to verify that \smash{$(\htheta(x), \, \hnu(x))$} be consistent. Here, we use the following assumption to guarantee consistency; this setup is general enough to cover both instrumental variables regression and quantile regression.
Finally, our consistency and Gaussianty results require using some specific settings for the trees from Algorithm (ref). In particular, we require that all trees be honest and regular in the sense of wager2015estimation, as follows. In order to satisfy the minimum split probability condition below, our implementation relies on the device of denil2014narrowing, whereby the number splitting variables considered at each step of the algorithm is random; specifically, we try $\min\cb{\max\cb{\text{Poisson$(m)$}, \, 1}, \, p}$ variables at each step, where $m > 0$ is a tuning parameter.
For generality, we set up Assumptions (ref)--(ref) in an abstract way. We end this section by showing that, in the context of our main problems of interest requiring Assumptions (ref)--(ref) is not particularly stringent. Further examples that satisfy the above assumptions will be discussed in Sections (ref) and (ref).
Given these assumptions, we are now ready to provide an asymptotic characterization of generalzed random forests. In doing so, we note that existing asymptotic analyses of regression forests, including mentch2016quantifying, scornet2015consistency and wager2015estimation, were built around the fact that regression forests are averages of regression trees grown over sub-samples, and can thus be analyzed as $U$-statistics hoeffding1948class. Unlike regression forest predictions, however, the parameter estimates \smash{$\htheta(x)$} from generalized random forests are not averages of estimates made by different trees; instead, we obtain \smash{$\htheta(x)$} by solving a single weighted moment equation as in (ref). Thus, existing proof strategies do not apply in our setting.
We tackle this problem using the method of influence functions as described by hampel1974influence; in particular, we are motivated by the analysis of newey1994kernel. The core idea of these methods is to first derive a sharp, linearized approximation to the local estimator \smash{$\htheta(x)$}, and then to analyze the linear approximation instead. In our setup, the influence function heuristic motivates a natural approximation \smash{$\ttheta^*(x)$} to \smash{$\htheta(x)$} as follows. Let $\rho_i^*(x)$ denote the influence function of the $i$-th observation with respect to the true parameter value $\theta(x)$, \smash{$\rho_i^*(x) := -\xi^\top V(x)^{-1} \psi_{\theta(x), \, \nu(x)}(O_i)$}. These quantities are closely related to the pseudo-outcomes (ref) used in our gradient tree splitting rule; the main difference is that, here, the $\rho_i^*(x)$ depend on the unknown true parameter values at $x$ and are thus inaccessible in practice. We use the $*$-superscript to remind ourselves of this fact.
Then, given any set of forest weights $\alpha_i(x)$ used to define the generalized random forest estimate \smash{$\htheta(x)$} by solving (ref), we can also define a pseudo-forest
which we will use as an approximation for \smash{$\htheta(x)$}. We note that, formally, this pseudo-forest estimate \smash{$\ttheta^*(x)$} is equivalent to the output of an (infeasible) regression forest with weights $\alpha_i(x)$ and outcomes $\theta(x) + \rho_i^*(x)$.
The upshot of this approximation is that, unlike \smash{$\htheta(x)$}, the pseudo-forest \smash{$\ttheta^*(x)$} is a $U$-statistic. Because $\ttheta^*(x)$ is a linear function of the pseudo-outcomes $\rho_i^*(x)$, we can write it as an average of pseudo-tree predictions \smash{$\ttheta^*(x) = \frac{1}{B} \sum_{b= 1}^B \ttheta^*_b(x)$} with \smash{$\ttheta^*_b(x) = \sum_{i = 1}^n \alpha_{ib}(x) \p{\theta(x) + \rho_i^*(x)}$}. Then, because each individual pseudo-tree prediction \smash{$\ttheta_b^*(x)$} is trained on a size-$s$ subsample of the training data drawn without replacement (see Section (ref)), \smash{$\ttheta^*(x)$} is an infinite-order $U$-statistic whose order corresponds to the subsample size, and so the arguments of mentch2016quantifying or wager2015estimation can be used to study the averaged estimator \smash{$\ttheta^*(x)$} using results about $U$-statistics hoeffding1948class,efron1981jackknife.
Following this proof strategy, the key difficulty is in showing that our influence-based statistic \smash{$\ttheta^*(x)$} is in fact a good approximation for \smash{$\htheta(x)$}. To do so, we start by establishing consistency of \smash{$\htheta(x)$} for $\theta(x)$ given our assumptions; we note that this is the only point in the paper where we use the fact that $\psi$ is the negative gradient of a convex loss as in Assumption (ref).
Building on this consistency result, we obtain a coupling of the desired type in Lemma (ref), the main technical contribution of this paper. We note that separating the analysis of moment estimators into a local approximation argument that hinges on consistency and a separate result that establishes consistency is standard; see, e.g., Chapter 5.3 of van2000asymptotic. The remainder of our analysis assumes that trees are grown on subsamples of size $s$ scaling as $s = n^\beta$ for some $\beta_{\min} < \beta < 1$, with
where $\pi$ and $\omega$ are as in Specification (ref). This scaling guarantees that the errors of forests are variance-dominated.
Given this coupling result, it now remains to study the asymptotics of \smash{$\ttheta^*(x)$}. In doing so, we re-iterate that \smash{$\ttheta^*(x)$} is exactly the output of an infeasible regression forest trained on outcomes $\theta(x) + \rho_i^*(x)$. Thus, the results of wager2015estimation apply directly to this object, and can be used to establish its Gaussianity. That we cannot actually compute \smash{$\ttheta^*(x)$} does not hinder an application of their results. Pursuing this approach, we find that given (ref), \smash{$\ttheta^*(x)$} and \smash{$\htheta(x)$} are both asymptotically normal. By extending the same argument, we could also show that the nuisance parameter estimates \smash{$\hnu(x)$} are consistent and asymptotically normal; however, we caution that the tree splits are not necessarily targeted to expressing heterogeneity in \smash{$\nu(x)$}, and so the resulting \smash{$\hnu(x)$} may not be particularly accurate in finite samples.
Theorem (ref) can also be used for statistical inference about $\theta(x)$. Given any consistent estimator $\hsigma_n(x)/\sigma_n(x) \rightarrow_p 1$ of the noise scale of \smash{$\htheta_n(x)$}, Theorem (ref) can be paired with Slutsky's lemma to verify that \smash{$\limn \, \mathbb{E}[\theta(x) \in (\htheta_n(x) \pm \Phi^{-1}(1 -\alpha/2) \hsigma_n(x))] = \alpha$}. Thus, in order to build asymptotically valid confidence intervals for $\theta(x)$ centered on $\htheta(x)$, it suffices to derive an estimator for \smash{$\sigma_n(x)$}.
In order to do so, we again leverage coupling with our approximating pseudo-forest \smash{$\ttheta^*(x)$}. In particular, the proof of Theorem (ref) implies that $\text{Var}[\ttheta^*(x)]/\sigma^2_n(x) \rightarrow_p 1$, and so it again suffices to study \smash{$\ttheta^*(x)$}. Moreover, from the definition of \smash{$\ttheta^*(x)$}, we directly see that
where \smash{$H_{n}(x; \, \theta, \, \nu) = \Var{\sum_{i = 1}^n \alpha_i(x) \psi_{\theta, \, \nu}(O_i)}$}. Thus, we propose building Gaussian confidence intervals using
where \smash{$\hV_n(x)$} and \smash{$\hH_n(x)$} are consistent estimators for the quantities in (ref).
The first quanitity $V(x)$ is a problem specific curvature parameter, and is not directly linked to forest-based methods. It is the same quantity that is needed to estimate variance of classical local maximum likelihood methods following newey1994kernel; e.g., for the instrumental variables problem described in Section (ref),
while for quantile regression, $V(x) = f_x(\theta(x))$. In both cases, several different strategies are available for estimating this term. In the case of instrumental variables forests, we suggest estimating the entries of (ref) using (honest and regular) regression forests.
The more interesting term is the inner variance term $H_n(x; \, \theta(x), \, \nu(x))$. To study this quantity, we note that the forest score \smash{$\Psi(\theta(x), \, \nu(x)) = \sum_{i = 1}^n \alpha_i(x) \psi_{\theta(x), \, \nu(x)}(O_i)$} is again formally equivalent to the output of a regression forest with weights $\alpha_i(x)$, this time with effective outcomes $\psi_{\theta(x), \, \nu(x)}(O_i)$. A number of proposals have emerged for estimating the variance of a regression forest, including work by sexton2009standard, mentch2016quantifying and wager2014confidence; and, in principle, any of these methods could be adapted to estimate the variance of $\Psi$. The only difficulty is that $\Psi$ depends on the true parameter values $(\theta(x), \, \nu(x))$, and so cannot directly be accessed in practice. Here, we present results based on a variant of the bootstrap of little bags algorithm (or noisy bootstrap) proposed by sexton2009standard. As a side benefit, we also obtain the first consistency guarantees for this method for any type of forest, including regression forests.
To motivate the bootstrap of little bags, we first note that building confidence intervals via half-sampling---whereby we evaluate an estimator on random halves of the training data to estimate its sampling error---is closely related to the bootstrap efron1982jackknife (throughout this section, we assume that $s \leq \lfloor n/2 \rfloor$). In our context, the ideal half-sampling estimator would be \smash{$\hH_n^{HS}(x)$} defined as
where \smash{$\Psi_{\hh}$} denotes a version of $\Psi$ computed only using all the possible trees that only rely on data from the half sample $\hh \subset \cb{1, \, ..., \, n}$ (specifically, in terms of Algorithm (ref), we only use trees whose full $\ii$-subsample is contained in $\hh$). If we could evaluate \smash{$\hH_n^{HS}(x)$}, results from efron1982jackknife suggest that it would be a good variance estimator for $\Psi$, but doing so is effectively impossible computationally as it would require growing very many forests.
Following sexton2009standard, however, we can efficiently approximate \smash{$\hH_n^{HS}(x)$} at almost no computational cost if we are willing to slightly modify our subsampling scheme. To do so, let $\ell \geq 2$ denote a little bag size and assume, for simplicity, that $B$ is an integer multiple of it. Then, we grow our forest as follows: First draw $g = 1, \, ..., \, B/\ell$ random half-samples $\hh_g \subset \cb{1, \, ..., \, n}$ of size $\lfloor n/2 \rfloor$, and then generate the subsamples $\ii_b$ used to build the forest in Algorithm (ref) such that $\ii_b \subseteq \hh_{\lceil b/\ell \rceil}$ for each $b = 1, \, ..., \, B$. In other words, we now generate our forest using little bags of $\ell$ trees, where all the trees in a given bag only use data from the same half-sample. sexton2009standard discuss optimal choices of $\ell$ for minimizing Monte Carlo error, and show that they depend on the ratio of the sampling variance of a single tree to that of the full forest.
The upshot of this construction is that we can now identify \smash{$\hH_n^{HS}(x)$} using a simple variance decomposition. Writing $\Psi_b$ for a version of $\Psi$ computed only using the $b$-th tree, we can verify that \smash{$\hH_n^{HS}(x)$} can be expressed in terms of the “between groups” and “within group” variance terms,
where $\mathbb{E}_{ss}$ denotes expectations over the subsampling mechanism while holding the data fixed. We define our feasible boostrap of little bags variance estimator \smash{$\hH_n^{BLB}(x)$} via a version of the above ANOVA decomposition that uses empirical moments and note that, given a large enough number of trees $B$, this converges to the ideal half-sampling estimator.
The result below verifies that, under the conditions of Theorem (ref), the optimal half-sampling estimator \smash{$\hH_n^{HS}(x)$} with plug-in values for \smash{$(\htheta(x), \, \hnu(x))$} as in (ref) consistently estimates the sampling variance of $\Psi(\theta(x), \, \nu(x))$. We have already seen above that the computationally feasible estimator \smash{$\hH_n^{BLB}(x)$} will match \smash{$\hH_n^{HS}(x)$} whenever $B$ is large enough and so, given any consistent estimator \smash{$\hV_n(x)$} for $V(x)$, we find that the confidence intervals built using (ref) will be asymptotically valid.
One challenge with the empirical moment estimator based on the above is that, if $B$ is small, the variance estimates \smash{$\hH_n^{BLB}(x)$} may be negative. In our software, we avoid this problem by using a Bayesian analysis of variance following, e.g., gelman2014bayesian, with an improper uniform prior for \smash{$\hH_n^{HS}(x)$} over $[0, \, \infty)$. When $B$ is large enough, this distinction washes out.
Our first application of generalized random forests is to the classical problem of non-parametric quantile regression. This problem has also been considered in detail by meinshausen2006quantile, who proposed a consistent forest-based quantile regression algorithm; his method also fits into the paradigm of solving estimating equations (ref) using random forest weights (ref). However, unlike us, meinshausen2006quantile does not propose a splitting rule that is tailored to the quantile regression context, and instead builds his forests using plain CART regression splits. Thus, a comparison of our method with that of meinshausen2006quantile provides a perfect opportunity for evaluating the value of our proposed method for constructing forest-based weights $\alpha_i(x)$ that are specifically designed to express heterogeneity in conditional quantiles.
Recall that, in the language of estimating equations, the $q$-th quantile $\theta_q(x)$ of the distribution of $Y$ conditionally on $X = x$ is identified via (ref), using the moment function \smash{$\psi_{\theta}(Y_i) = q \1\p{\cb{Y_i > \theta}} - (1 - q) \1\p{\cb{Y_i \leq \theta}}$}. Plugging this moment function into our splitting scheme (ref) gives us pseudo-outcomes \smash{$\rho_i = \1(\{Y_i > \htheta_{q, P(X_i)}\})$}, where \smash{$\htheta_{q, P(X_i)}$} is the $q$-th quantile of the parent node $P(X_i)$ containing $X_i$, up to a scaling and re-centering that do not affect the subsequent regression split on these pseudo-outcomes. In other words, gradient-based quantile regression trees try to separate observations that fall above the $q$-th quantile of the parent from those below it.
We compare our method to that of meinshausen2006quantile in Figure (ref). In the left panel, we have a mean shift in the distribution of $Y_i$ conditional on $X_i$ at $(X_i)_1 = 0$, and both methods are able to pick it up as expected. However, in the right panel, the mean of $Y$ given $X$ is constant, but there is a scale shift at $(X_i)_1 = 0$. Here, our method still performs well, as our splitting rule targets changes in the quantiles of the $Y$-distribution. However, the method of meinshausen2006quantile breaks down completely, as it relies on CART regression splits that are only sensitive to changes in the conditional mean of $Y$ given $X$. We also note that generalized random forests produce somewhat smoother sample paths than the method of meinshausen2006quantile; this is due to our use of honesty as described in Section (ref). If we run generalized random forests without honesty, then our method still correctly identifies the jumps at $x = 0$, but has sample paths that oscillate locally just as much as the baseline method. The purpose of this example is not to claim that our variant of quantile regression forests built using gradient trees is always superior to the method of meinshausen2006quantile that uses regression-based splitting to obtain the weights $\alpha_i(x)$; rather, we found that, our splitting rule is specifically sensitive to quantile shifts in a way that regression splits are not---and, moreover, deriving our splitting rule was fully automatic given the generalized random forest formalism.
In several applications, we want to estimate multiple quantiles at the same time. For example, in Figure (ref), we estimate at $q = 0.1, \, 0.5, \, 0.9$. Estimating different forests for each quantile separately would be undesirable for many reasons: It would be computationally expensive, and there is a risk that quantile estimates might cross in finite samples due to statistical noise. Thus, we need to build a forest using a splitting scheme that is sensitive to changes at any of our quantiles of interests. Here, we use a simple heuristic inspired by our relabeling transformation. Given a set of quantiles of interest $q_1 < ... < q_k$, we first evaluate all these quantiles \smash{$\htheta_{q_1, P(X_i)} \leq ... \leq \htheta_{q_k, P(X_i)}$} in the parent node, and label $i$-th point by the interval \smash{$[\htheta_{q_{j - 1}, P(X_i)}, \, \htheta_{q_{j}, P(X_i)})$} it falls into. Then, we choose the split point using a multiclass classification rule that classifies each observation into one of the intervals.
Next, we consider conditional average partial effect estimation under exogeneity; procedurally, the statistical task is equivalent to solving linear regression problems conditionally on features. Suppose that we observe samples $(X_i, \, Y_i, \, W_i) \in \xx \times \RR \times \RR^q$, and posit a random effects model \smash{$Y_i = W_i \cdot b_i + \varepsilon_i$}, \smash{$\beta(x) = \EE{b_i \cond X_i = x}$}. Our goal is to estimate $\theta(x) = \xi \cdot \beta(x)$ for some contrast $\xi \in \RR^p$. If $W_i \in \cb{0, \, 1}$ is a treatment assignment, then $\beta(x)$ corresponds to the conditional average treatment effect.
In order for the average effect $\beta(x)$ to be identified, we need to make certain distributional assumptions. Here, we assume that the $W_i$ are exogenous, i.e., independent of the unobservables conditionally on $X_i$: \smash{$\cb{b_i, \, \varepsilon_i} \indep W_i \cond X_i$}. If $W_i$ is a binary treatment, this condition is equivalent to the unconfoundedness assumption used to motivate propensity score methods rosenbaum1983central. When exogeneity does not hold, more sophisticated identification strategies are needed (see following section).
Our parameter of interest $\theta(x) = \xi \cdot \beta(x)$ is identified by (ref) with \smash{$\psi_{\beta(x), \, c(x)}(Y_i, \, W_i) = (Y_i - \beta(x) \cdot W_i - c(x)) (1 \ \ W_i^\top)^\top$} where $c(x)$ is an intercept term; this can also be written more explicitly as \smash{$\theta(x) = \xi^\top \Var{W_i \cond X_i = x}^{-1} \Cov{W_i, \, Y_i \cond X_i = x}$}. Given forest weights $\alpha_i(x)$ as in (ref), the induced estimator \smash{$\htheta(x)$} for $\theta(x)$ is
where \smash{$\bW_\alpha = \sum \alpha_i(x) W_i$} and \smash{$\bY_\alpha = \sum \alpha_i(x) Y_i$}, and we write \smash{$v^{\otimes 2} = vv^\top$}.
Generalized random forests provide us with a quasi-automatic framework for getting the weights \smash{$\alpha_i(x)$} needed in (ref); all that needs to be done is to compute the pseudo-outcomes \smash{$\rho_i$} from (ref) used for recursive partitioning. We use (ref) and, for every parent $P$ and each observation $i$ with $X_i \in P$ set
where now \smash{$\bW_P$} and \smash{$\bY_P$} stand for averages taken over the parent $P$, and \smash{$\hbeta_P$} is the least-squares regression solution of $Y_i$ on $W_i$ in the parent. Note that the matrix inverse \smash{$A_P^{-1}$} only needs to be evaluated once per parent node.
Checking the conditions required in Section (ref), note that Assumption (ref) holds whenever the functions $\EE{Y_i \cond X_i = x}$, $\EE{W_i \cond X_i = x}$, $\Cov{Y_i, \, W_i \cond X_i = x}$ and $\Var{W_i \cond X_i = x}$ are all Lipschitz in $x$, Assumption (ref) holds provided that $\Var{W_i \cond X_i = x}$ in invertible, while Assumptions (ref)--(ref) hold by construction. Thus, Theorem (ref) in fact applies in this setting.
The above construction allows for asymptotically valid inference for $\theta(x)$, but the performance of the forests can in practice be improved by first regressing out the effect of the features $X_i$ on all the outcomes separately. Writing \smash{$y(x) = \mathbb{E}[Y_i \cond X = x]$} and \smash{$w(x) = \mathbb{E}[W_i \cond X = x]$} for the conditional marginal expectations of $Y_i$ and $W_i$ respectively, define centered outcomes \smash{$\tY_i = Y_i - \hy^{(-i)}\p{X_i}$} and \smash{$\tW_i = W_i - \hw^{(-i)}\p{X_i}$}, where \smash{$\hy^{(-1)}\p{X_i}$}, etc., are leave-one-out estimates of the marginal expectations, computed without using the $i$-th observation. We then run a forest using centered outcomes \smash{$\{\tY_i, \, \tW_i\}_{i = 1}^n$} instead of the original \smash{$\{Y_i, \, W_i\}_{i = 1}^n$}.
In order to justify this transformation, we note if there is any set $\set \subseteq \xx$ over which $\beta(x)$ is constant (and so $\theta(x)$ is also constant), the following expression also identifies $\theta(x)$ for any $x \in \set$:
Thus, if we locally center the $Y_i$ and the $W_i$ before running our forest, the estimator (ref) has the potential to be robust to confounding effects even when the weights $\alpha_i(x)$ are not sharply concentrated around $x$. Similar orthogonalization ideas have proven to be useful in many statistical contexts chernozhukov2016double,newey1994asymptotic,neyman1979c; in particular, robinson1988root showed that if we have access to a neighborhood $\set$ over which $\beta(x) = \beta_{\set}$ is constant, then the moment condition (ref) induces a semiparametrically efficient estimator for $\theta_{\set} = \xi \cdot \beta_{\set}$.
We note that if we ran a forest with any deterministic centering scheme, i.e., we used \smash{$\tY_i = Y_i - \hy(X_i)$} for any Lipschitz function $\hy(X_i)$ that does not depend on the data, etc., then the theory developed in Section (ref) would allow for valid inference about $\theta(x)$ (in particular, we do not need to assume consistency of $\hy(X_i)$). Moreover, we could also emulate this result by using a form of $k$-fold cross-fitting chernozhukov2016double,schick1986asymptotically. In the context of forests, it is much more practical to carry out residualization via leave-one-out prediction than via $k$-fold cross-fitting, because leave-one-out prediction in forests is computationally cheap breiman2001random; however, a practitioner wanting to use results that are precisely covered by theory may prefer to use cross-fitting for centering.
When $W_i \in \cb{0, \, 1}$ is a binary treatment assignment, the present setup is equivalent to the standard problem of heterogeneous treatment effect estimation under unconfoundedness. Heterogeneous treatment effect estimation via tree-based methods has received considerable attention in the recent literature: athey2016recursive and su2009subgroup develop tree-based methods for subgroup analysis, hill2011bayesian studies treatment effect estimation via Bayesian additive regression trees chipman2010bart, and wager2015estimation propose a causal forest procedure that is very nearly a special case of our generalized random forests. The main interest of our method is in how it can handle situations for which no comparable methods exist, such as instrumental variables regression as discussed below. Here, however, we briefly discuss how some concepts developed as a part of our more general approach directly improve the performance of causal forests.
The closest method to ours is Procedure 1 of wager2015estimation, which is almost equivalent to a generalized random forest without centering, the only substantive differences being that they split using the exact loss criterion (ref) rather than our gradient-based loss criterion (ref), and let each tree compute its own treatment effect estimate rather than using the weighting scheme from Section (ref) (these methods are exactly equivalent for regression forests, but not for causal forests). wager2015estimation also consider a second approach, Procedure 2, that obtains its neighborhood function via a classification forest on the treatment assignments $W_i$.
A weakness of the methods in wager2015estimation, as they note in their discussion, is that these two procedures have different strengths---Procedure 1 is more sensitive to changes in the treatment effect function, while Procedure 2 is more robust to confounding---but the hard coded nature of these methods made it difficult to reconcile their relative advantages. Conversely, given the framing of generalized random forests via estimating equations, it is “obvious” that we can leverage best practices from the literature on estimating equations and orthogonalize our moment conditions by regressing out the main effect of $X_i$ on $W_i$ and $Y_i$ as in robinson1988root.
To illustrate the value of orthogonalization, we revisit a simulation of wager2015estimation where $X_i \sim U([0, \, 1]^p)$, $W_i \cond X_i \sim \text{Bernoulli}(e(X_i))$, and $Y_i \cond X_i, \, W_i \sim \nn\p{m(X_i) + (W_i - 0.5) \tau(X_i), \, 1}$. The authors consider two different simulation settings: One with no confounding, $m(x) = 0$ and $e(x) = 0.5$, but with treatment heterogeneity \smash{$\tau\p{x} = \varsigma\p{x_1} \, \varsigma\p{x_2}$}, \smash{$\varsigma\p{u} = 1 + {1}\,/\,({1 + e^{-20\p{u - 1/3}}})$}, and second with no treatment effect, $\tau(x) = 0$, but with confounding, \smash{$e(x) = \frac{1}{4} \p{1 + \beta_{2, \, 4}(x_3)}, \ \ m(x) = 2x_3 - 1$}, where $\beta_{a, \, b}$ is the $\beta$-density with shape parameters $a$ and $b$. We also consider a third setting with both heterogeneity and confounding, that combines $\tau(\cdot)$ from the first setting with $m(\cdot)$ and $e(\cdot)$ from the second. For the first setting, wager2015estimation used their Procedure 1, whereas for the second they used Procedure 2, while noting that it is unfortunate that the practitioner is forced to choose one procedure or the other.
Results presented in Table (ref) are reassuring, suggesting that generalized random forests with centering do well under both settings, and can better handle the case with both confounding and treatment heterogeneity than either of the other two procedures. In contrast, Procedure 1 of Wager and Athey does poorly with pure confounding, whereas Procedure 2 of Wager and Athey is good in the pure confounding setting, but does poorly with strong heterogeneity; this is as expected, noting the design of both methods.
In many applications, we want to measure the causal effect of an intervention on an outcome, all while recognizing that the intervention and the outcome may also be tied together through non-causal pathways, thus ruling out the exogeneity assumption made above. One approach in this situation is to rely on instrumental variables (IV) regression, where we find an auxiliary source of randomness that can be used to identify causal effects.
For example, suppose we want to measure the causal effect of child rearing on a mother's labor-force participation. It is well known that, in the United States, mothers with more children are less likely to work. But how much of this link is causal, i.e., some mothers work less because they are busy raising children, and how much of it is merely due to confounding factors, e.g., some mothers have preferences that both lead them to raise more children and be less likely to participate in the labor force? Understanding these effects may be helpful in predicting the value of programs like subsidized daycare that assist mothers' labor force participation while they have young children.
To study this question, angrist1998children found a source of auxiliary randomness that can be used to distinguish causal versus correlational effects: They found that, in the United States, parents who already have two children of mixed sexes, i.e., one boy and one girl, will have fewer kids in the future than parents whose first two children were of the same sex. Assuming that the sexes of the first two children in a family are effectively random, this observed preference for having children of both sexes provides an exogenous source of variation in family size that can be used to identify causal effects: If the mixed sex indicator is unrelated to the mother's propensity to work for a fixed number of children, then the effect of the mixed sex indicator on the observed propensity to work can be attributed to its effect on family size. The instrumental variable estimator normalizes this effect by the effect of mixed sex on family size, so that the normalized estimate is a consistent estimate of the treatment effect of family size on work. Other classical uses of instrumental variables regression include measuring the impact of military service on lifetime income by using the Vietnam draft lottery as an instrument angrist1990lifetime, and measuring the extent to which 401(k) savings programs crowd out other savings, using eligibility for 401(k) savings programs as an instrument abadie2003semiparametric,poterba1996retirement.
Classical instrumental variables regression seeks a global understanding of the treatment effect, e.g., on average over the whole U.S. population, does having more children reduce the labor force participation of women? Here, we use forests to estimate heterogeneous treatment effects: We might ask how the causal effect of child rearing varies with a mother's age and socioeconomic status.
We observe $i = 1, \, ..., \, n$ independent and identically distributed subjects, each of whom has features $X_i \in \xx$, an outcome $Y_i \in \RR$, a received treatment $W_i \in \cb{0, \, 1}$, and an instrument $Z_i \in \cb{0, \, 1}$. We believe that the outcomes $Y_i$ and received treatment $W_i$ are related via a structural model \smash{$Y_i = \mu\p{X_i} +\tau\p{X_i} W_i + \varepsilon_i$}, where $\tau(X_i)$ is understood to be the causal effect of $W_i$ on $Y_i$, and $\varepsilon_i$ is a noise term that may be positively correlated with $W_i$. Because $\varepsilon_i$ is correlated with $W_i$, standard regression analyses will not in general be consistent for $\tau(X_i)$, and we need to use the instrument $Z_i$. If $Z_i$ is independent of $\varepsilon_i$ conditionally on $X_i$ then, provided that $Z_i$ has an influence on the received treatment $W_i$, i.e., that the covariance of $Z_i$ and $W_i$ conditionally on $X_i = x$ is non-zero, the treatment effect $\tau(x)$ is identified via \smash{$\tau(x) = \Cov{Y_i, \, Z_i \cond X_i = x} \,\big/\, \Cov{W_i, \, Z_i \cond X_i = x}$}. We can then estimate $\tau(x)$ by via moment functions \smash{$ \EE{Z_i\p{Y_i - W_i \, \tau(x) - \mu(x)} \cond X_i = x} = 0$} and \smash{$\EE{Y_i - W_i \, \tau(x) - \mu(x) \cond X_i = x} = 0$}, where the intercept $\mu(x)$ is a nuisance parameter. If we are not willing to assume that every individual $i$ with features $X_i = x$ has the same treatment effect $\tau(x)$, then heterogeneous instrumental variables regression allows us to estimate a (conditional) local average treatment effect abadie2003semiparametric,imbens1994identification.
We then use our formalism to derive a forest that is targeted towards estimating causal effects identified via conditional two-stage least squares. Gradient-based labeling (ref) yields pseudo-outcomes for every parent node $P$ and each observation $i$ with $X_i \in P$, \smash{$\rho_i = \p{Z_i - \bZ_P} \p{\p{Y_i - \bY_P} - \p{W_i - \bW_P} \htau_P}$}, where $\bY_P, \, \bW_P, \, \bZ_P$ are moments in the parent node, and $\htau_P$ is a solution to the estimating equation in the parent. Given these pseudo-outcomes, the tree executes a CART regression split on the $\rho_i$ as usual. Finally, we obtain personalized treatment effect estimates $\htau(x)$ by solving (ref) with forest weights (ref).
To verify that Theorem (ref) holds in this setting, we note that Assumption (ref) holds whenever the conditional moment functions $\EE{W_i \cond X_i = x}$, $\EE{Y_i \cond X_i = x}$, $\EE{Z_i \cond X_i = x}$, $\EE{W_iZ_i \cond X_i = x}$ and $\EE{Y_iZ_i \cond X_i = x}$ are all Lipschitz continuous in $x$, while Assumption (ref) holds whenever the instrument is correlated with received treatment (i.e., the instrument is valid). Assumptions (ref)--(ref) hold thanks to the definition of $\psi$.
As in Section (ref), we center our procedure using the transformation of robinson1988root, and regress out the marginal effects of $X_i$ first. Writing \smash{$y(x) = \mathbb{E}[Y_i \cond X = x]$}, \smash{$w(x) = \mathbb{E}[W_i \cond X = x]$}, and \smash{$z(x) = \mathbb{E}[Z_i \cond X = x]$}, we compute conditionally centered outcomes by leave-one-out estimation \smash{$\tY_i = Y_i - \hy^{(-i)}\p{X_i}$}, \smash{$\tW_i = W_i - \hw^{(-i)}\p{X_i}$} and \smash{$\tZ_i = Z_i - \hz^{(-i)}\p{X_i}$}, \sloppy{and then run the full instrumental variables forest using centered outcomes \smash{$\{\tY_i, \, \tW_i, \, \tZ_i\}_{i = 1}^n$}. We recommend working with centered outcomes by default, and we do so in our simulations. Our package grf provides the option of making this transformation automatically, where \smash{$\hy^{(-i)}\p{X_i}$}, \smash{$\hw^{(-i)}\p{X_i}$} and \smash{$\hz^{(-i)}\p{X_i}$} are first estimated using 3 separate regression forests.}
There is a rich literature on non-parametric instrumental variables regression. The above approach generalizes classical approaches based on kernels or series estimation abadie2003semiparametric,su2013local,wooldridge2010econometric. In other threads, darolles2011nonparametric and newey2003instrumental study instrumental variables models that generalize the conditionally linear treatment model and allowing for non-linear effects, and hartford2016counterfactual develop deep learning tools. belloni2012sparse consider working with high-dimensional instruments $Z_i$.
Appendix (ref) has a simulation study for IV forests, comparing them to nearest-neighbor and series regression. We find our method to perform well relative to these baselines, and centering to be helpful. We also evaluate coverage of the bootstrap of little bag confidence intervals.
We now revisit our motivating example discussed at the beginning of Section (ref). We follow angrist1998children in constructing our dataset, and study a sample of $n = 334,535$ married mothers with at least 2 children (1980 census data), based on the following quantities: The outcome $Y_i$ is whether the mother did not work in the year preceding the census, the received treatment $W_i$ is whether the mother had 3 or more children at census time, and the instrument $Z_i$ measures whether or not the mother's first two children were of different sexes. Based on this data, angrist1998children estimated the local average treatment effect of having a third child among mothers with at least two children. In our sample, \smash{$\hCov{W, \, Z} = 1.6 \cdot 10^{-2}$}, while \smash{$\hCov{Y, \, Z} = 2.1 \cdot 10^{-3}$}, leading to a 95% confidence interval for the local average treatment effect $\tau \in (0.14 \pm 0.054)$ using the R function ivreg kleiber2008applied. Thus, it appears that having a third child reduces women's labor force participation on average in the US. angrist1998children conduct extensive sensitivity analysis to corroborate the plausibility of this identification strategy.
We seek to extend this analysis by fitting heterogeneity on several covariates, including the mother's age at the birth of her first child, her age at census time, her years of education and her race (black, hispanic, other), as well as the father's income. Formally, our analysis identifies a conditional local average treatment effect $\tau(x)$ abadie2003semiparametric,imbens1994identification.
Results from a generalized random forest analysis are presented in Figure (ref). These results suggest that the observed treatment effect is driven by mothers whose husbands have a lower income. Such an effect would be intuitively easy to justify: it seems plausible that mothers with wealthier husbands can afford to hire help in raising their children, and so can choose whether or not to work based on other considerations. That being said, we caution that the father's income was measured in the census, so there is potentially an endogeneity problem: perhaps a mother's choice not to work after having a third child enables the husband to earn more. Ideally, we would have wanted to measure the husband's income at the time of the second child's birth, but we do not have access to this measurement in the present data. Moreover, the confidence intervals in Figure (ref) are rather wide, attesting to the importance of formal asymptotic theory when using forest-based methods for instrumental variables regression.
We introduced generalized random forests as a versatile method for adaptive, local estimation in a wide variety of statistical models. We discussed our method in the contexts of quantile regression and heterogeneous treatment effect estimation, and our approach also applies to a wide variety of other settings, such as demand estimation or panel data analysis. Our software, grf, is implemented in a modular way that should enable users to implement splitting rules motivated by new statistical questions.
Many of the remaining challenges with generalized random forests are closely related to those with standard nonparametric methods for local likelihood estimation. In particular, as discussed above, our confidence interval construction relies on undersmoothing to get valid asymptotic coverage (without undersmoothing, the confidence intervals account for sampling variability of the forest, but do not capture bias). Developing a principled way to bias-correct our confidence intervals, and thus avoid the need for undersmoothing, would be of considerable interest both conceptually and in practice. Moreover, again like standard methods, forests can exhibit edge effects whereby the slope of our estimates \smash{$\htheta(x)$} may taper off as we approach the edge of $\xx$-space, even when the true function $\theta(x)$ keeps changing. Finding an elegant way to deal with such edge effects could improve the quality of the confidence intervals provided by generalized random forests.
We are grateful to Jerry Friedman for first recommending we take a closer look at splitting rules for quantile regression forests, to Will Fithian for drawing our attention to connections between our early ideas and gradient boosting, to Guido Imbens for suggesting the local centering scheme in Section (ref), to the associate editor and three anonymous referees for helpful suggestions, and to seminar participants at the Atlantic Causal Inference Conference, the BIRS Workshop on the Interface of Machine Learning and Statistical Inference, the California Econometrics Conference, Ca'Voscari University of Venice, Columbia, Cornell, the Econometric Society Winter Meetings, EPFL, the European University Institute, INFORMS, Kellogg, the Microsoft Conference on Digital Economics, the MIT Conference on Digital Experimentation, Northwestern, Toulouse, Triangle Computer Science Distinguished Lecture Series, University of Chicago, University of Illinois Urbana--Champaign, University of Lausanne, and the USC Dornsife Conference on Big Data in Economics.