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.
95,648 characters · 20 sections · 35 citation commands
Balancing Flexibility and Interpretability: A Conditional Linear Model Estimation via Random Forest
\footnotetext[1]{ Department of Statistics, University of California, Davis. } \footnotetext[2]{ Department of Economics, University of Illinois, Urbana-Champaign. } \let\arabic{footnote}\relax \footnotetext[1]{ \textsuperscript{*}Corresponding author: [email removed]} }
\setcounter{page}{0}\thispagestyle{empty}
\def\spacingset#1{ {#1}} \spacingset{1}
\spacingset{1.5}
In economics and related social sciences, epidemiology and medicine, psychology, and many other areas, estimating the partial effects (causal or not) of one or more factors on a target variable is extremely important. In general, partial effects estimation is conducted either on parametric models with strong functional-form assumptions or in overly simplified semi-parametric alternatives. These models provide interpretive clarity and computational efficiency but often impose restrictive assumptions that may inadequately capture the complexity of real-world data. In recent years, machine learning (ML) methods have broadened the statistician/econometrician’s toolkit by offering more flexible approaches to modeling relationships within potentially high-dimensional datasets without imposing stringent parametric assumptions; see, for example, sAgI2019 and rMeMmM2023 for recent review papers. Despite these advances, the challenge of balancing flexibility with interpretability remains a significant issue for applied researchers.
This paper presents a locally linear model that addresses the challenge of integrating the flexibility of machine learning (ML) methods with the interpretability of linear models. The proposed model aligns with several well-established parametric and semi-parametric specifications in the literature, including switching regression mgD1969,smGreQ1972, varying-coefficient models, and additive models tHrT1993,rCrsT1993b. It also incorporates recent advancements in machine learning, such as Random Forests (RF) introduced by breiman2001random, Generalized Random Forests (GRF) presented by sAjTsW2019, and local linear forests (LLF) discussed in rFjTsAsW2021. Our central contribution is a robust framework that enables the nonparametric estimation of heterogeneous partial effects of one or more variables of interest on an outcome. Our approach is intuitive and computationally simple, achieved through an adaptation of the RF method.
Let $Y$ be a response random variable and define $h(\boldsymbol{x},\boldsymbol{z}):=\E(Y|\boldsymbol{X}=\boldsymbol{x},\boldsymbol{Z}=\boldsymbol{z})$, where $h(\cdot)$ is an unknown Lipschitz continuous function of two vectors of random covariates $\boldsymbol{X}$ and $\boldsymbol{Z}$. Suppose the goal is to estimate $h(\boldsymbol{x},\boldsymbol{z})$ for arbitrary data points $\boldsymbol{x}$ and $\boldsymbol{z}$, as well as the partial effects $\partial h(\boldsymbol{x},\boldsymbol{z})/\partial \boldsymbol{x}$ assuming $h$ is differentiable with respect to its first argument. Modern machine learning methods provide a range of nonparametric algorithms to estimate $h(\boldsymbol{x},\boldsymbol{z})$, enabling greater flexibility in modeling intricate relationships. However, these models often lack interpretability, and the computation of partial effects can be both intensive and non-trivial, particularly in the context of deep learning models.\footnote{Deep learning methods usually require the definition of too many hyperparameters, the network architecture, and also require large amounts of data for effective training.} Moreover, conducting inference on partial effects within the framework of general machine learning models remains an open problem, necessitating further research and development to enhance methodological rigor and applicability in statistical analysis.
This paper presents an alternative approach rooted in contemporary machine learning literature, specifically aimed at enhancing interpretability. We consider covariates $\boldsymbol{X}$ and $\boldsymbol{Z}$, which are not necessarily mutually independent, with a focus on the estimation $\partial h(\boldsymbol{x},\boldsymbol{z})/\partial \boldsymbol{x}$. We propose a model in which the conditional expectation of $Y$, given $\boldsymbol{X}$ and $\boldsymbol{Z}$, is represented as follows: \[ h(\boldsymbol{x},\boldsymbol{z}):=\E(Y|\boldsymbol{X}=\boldsymbol{x},\boldsymbol{Z}=\boldsymbol{z})=\boldsymbol\beta(\boldsymbol{z})^{\ensuremath{\mathsf{T}}}\boldsymbol{x}. \]
This is a varying-coefficient model where the coefficients of the linear relationship between $\boldsymbol{X}$ and $Y$ vary as a function of $\boldsymbol{Z}$. Importantly, the partial effect of $\boldsymbol{X}$ on $Y$, $\partial h(\boldsymbol{x},\boldsymbol{z})/\partial \boldsymbol{x}$, is given directly by $\boldsymbol\beta(\boldsymbol{z})$, and thus exhibits heterogeneity across different values of $\boldsymbol{Z}$. Estimating the map $\boldsymbol{z}\mapsto \boldsymbol\beta(\boldsymbol{z})$ is equivalent to recovering the heterogeneous partial effects, providing a natural and intuitive model interpretation. Clearly, this is only the case when $\boldsymbol{Z}$ and $\boldsymbol{X}$ do not share the same covariates. We consider the estimation of $\boldsymbol{\beta}(\boldsymbol{z})$ by modification of the Random Forest (RF) method, where the trees have a linear model on each of their leaves. If $\boldsymbol{X}$ is a fixed scalar, the model equals the Random Forest (RF) method, where the trees have an intercept model on each of their leaves.
This model presents two primary advantages. First, it facilitates flexible, data-driven estimation of complex relationships, all while ensuring the interpretability of partial effects, which is often a pivotal focus in applied research. Second, by integrating the flexibility of the Random Forest (RF) algorithm within a locally linear framework, we effectively capture heterogeneity in partial effects in a computationally feasible manner. This characteristic renders the approach particularly well-suited for large-scale empirical applications.
In many empirical applications, the dimensionality of $\boldsymbol{X}$ is expected to be small relative to the sample size (it is not uncommon to have an univariate $\boldsymbol{X}$), simplifying the application of our model. Partial effects of covariates can be estimated directly using modern machine learning and nonparametric techniques. We choose the RF algorithm for its robustness and flexibility, mainly because it requires fewer hyperparameter choices, mitigating the risks of overfitting and cherry-picking. Additionally, RF's efficient estimation algorithms make it computationally attractive. While alternative approaches like deep learning offer comparable flexibility, they are more sensitive to hyperparameter tuning and initial conditions, require larger datasets, and are computationally intensive. By contrast, RF offers a more stable and efficient solution.
The concept of introducing nonlinearity through varying coefficients in linear models is well-established. It originated in threshold regression models by mgD1969, rQ1972, smGreQ1972, nmK1978, hTksL1980, and later work by rsT1989. These studies focus on sharp parameter changes based on a univariate threshold variable, capturing structural shifts in the data. Smooth transitions (shifts) were introduced by ksChT1986a and tT1994a, but these models relied on a single transition variable. Extensions to multiple transition variables inspired by the neural network (NN) literature were explored by mcMaV2000, mcMaV2005, and sMcPmcM2004, though these approaches were parametric and limited to low-dimensional settings. Nonparametric alternatives have been proposed by tHrT1993, jFwZ1999, and zCjFqY2000. However, estimating these models becomes challenging with an increase in the number of covariates. In contrast, our semi-parametric approach based on random forests accommodates a large set of covariates, provided the dimensionality remains smaller than the sample size.
Estimating regression trees with linear models in the terminal nodes is not new and is nested in the more general Generalized Random Forest model by sAjTsW2019. rFjTsAsW2021 formalized the locally linear random forests and recommended a local Ridge regression to estimate the parameters. Nevertheless, this paper differs and complements the ones cited above in a few ways. First, we estimate the local linear models with the usual ordinary least-squares method, making our algorithm simpler than the ones in sAjTsW2019 and rFjTsAsW2021. Second, as our primary focus is on the estimation of the partial effects $\boldsymbol{\beta}(\boldsymbol{z})$ and not $\E(Y|\boldsymbol{X}=\boldsymbol{x},\boldsymbol{Z}=\boldsymbol{z})$, we complement the previous papers by deriving a new consistency and asymptotic normality results for $\widehat{\boldsymbol\beta}(\boldsymbol{z})$. We also worked out the rates of convergence and derived the asymptotic covariance matrix for the estimator of $\widehat{\boldsymbol\beta}(\boldsymbol{Z})$. Third, we proposed a consistent estimator for the covariance matrix of $\widehat{\boldsymbol\beta}(\boldsymbol{Z})$. Fourth, unlike most papers in the literature, we show that our results are also valid when $\boldsymbol{Z}$ contains discrete random variables. Note that having discrete-valued elements of $\boldsymbol{Z}$ is not the same as introducing interaction effects between $\boldsymbol{X}$ and dummy variables constructed from the classes in $\boldsymbol{Z}$. Our approach allows the partial effect heterogeneity to be determined by unknown interactions among the elements of $\boldsymbol{Z}$. Finally, based on our new convergence results, we derived two tests to conduct inference on the partial effects and test whether the partial effects are homogeneous.
Following this introduction, the paper is structured as follows. In Section (ref), we define the model, outline assumptions related to the data-generating mechanism, and present special cases of our proposal. In Section (ref), we present the theoretical results of this paper. More specifically, Section (ref) details the main convergence results, while Sections (ref) and (ref) provide descriptions of the specification tests proposed in the paper. The case of discrete random variables is examined in Section (ref). Monte Carlo simulations are illustrated in Section (ref), and empirical samples are presented in Section (ref). Finally, Section (ref) wraps up the paper. All technical derivations are provided in the Supplemental Material.
Vectors are denoted by bold lowercase $\boldsymbol{x}$ and matrices by bold uppercase $\boldsymbol{M}$. $\boldsymbol{x}^{\ensuremath{\mathsf{T}}}$ denotes the transpose of the vector $\boldsymbol{x}$. Similarly $\boldsymbol{M}^{\ensuremath{\mathsf{T}}}$ denotes the transpose for matrix $\boldsymbol{M}$. We write $\|\boldsymbol{x}\|_p$ for $p\in[1,\infty]$ to denote the $\ell^p$-norm if $\boldsymbol{x}$ is a (possibly random) vector or the induced operator $\ell^p$--$\ell^p$-norm if $\boldsymbol{M}$ is a matrix. For a matrix $\boldsymbol{M}$, we write $\|\boldsymbol{M}\|_{\max}$ for the maximum absolute entry and $\|\boldsymbol{M}\|_\ensuremath{\mathrm{F}}$ for the Frobenius norm. We denote positive semi-definiteness by $\boldsymbol{M} \succeq 0$ and write $\boldsymbol{I}_d$ for the $d \times d$ identity matrix.
For scalar sequences $x_n$ and $y_n$, we write $x_n \lesssim y_n$ if there exists a positive constant $C$ such that $|x_n| \leq C |y_n|$ for sufficiently large $n$. We write $x_n \asymp y_n$ to indicate both $x_n \lesssim y_n$ and $y_n \lesssim x_n$. Similarly, for random variables $X_n$ and $Y_n$, we write $X_n \lesssim_\ensuremath{\mathbb{P}} Y_n$ if for every $\varepsilon > 0$ there exists a positive constant $C$ such that $\ensuremath{\mathbb{P}}(|X_n| \geq C |Y_n|) \leq \varepsilon$, and write $X_n \overset{\ensuremath{\mathbb{P}}}{\longrightarrow} X$ and $X_n \overset{d}{\longrightarrow} X$ for limits in probability and in distribution, respectively. For real numbers $a$ and $b$ we use $a \land b = \max\{a,b\}$ and $a \lor b = \min\{a,b\}$.
We consider the following model.
The local-linear model (ref) has the advantage of having the marginal treatment effect built-in since $h(\boldsymbol{x},\boldsymbol{z})/\partial \boldsymbol{x}=\boldsymbol\beta(\boldsymbol{z})$. At the same time, the model is flexible enough in terms of the covariates $\boldsymbol{Z}$ to accommodate complex heterogeneous partial effects, which include higher-order interactions between $\boldsymbol{X}$ and $\boldsymbol{Z}$. Model (ref) nests several notable cases of interest.
Since $\E[Y|\boldsymbol{X},\boldsymbol{Z}]$ is the minimizer of $\E[(Y-f(\boldsymbol{X},\boldsymbol{Z}))^2]$ over the class of functions $f$ such that $\E[f(\boldsymbol{X},\boldsymbol{Z})^2]<\infty$, we can explicitly characterize the function $\boldsymbol{z}\mapsto \beta(\boldsymbol{z})$ such that $\boldsymbol\beta(\boldsymbol{z}) = \boldsymbol{\Omega}(\boldsymbol{z})^{-1}\boldsymbol{\gamma}(\boldsymbol{z})$, where $\boldsymbol{\Omega}(\boldsymbol{z}):= \E[\boldsymbol{X}\boldsymbol{X}^{\ensuremath{\mathsf{T}}}|\boldsymbol{Z}=\boldsymbol{z}]$ and $\boldsymbol{\gamma}(\boldsymbol{z}):=\E[\boldsymbol{X}Y|\boldsymbol{Z}=\boldsymbol{z}]$, and provided that $\boldsymbol{\Omega}(\boldsymbol{z})$ is almost sure positive definite.
Given a random sample $\{(Y_i,\boldsymbol{X}_i,\boldsymbol{Z}_i):1\leq i\leq n\}$ of $(Y,\boldsymbol{X},\boldsymbol{Z})$ we propose to estimate (ref) by a modification of the Random Forest procedure (refer to Algorithm (ref)). We start by defining for any nonempty set of indices $\mathcal{I}\subseteq [n]$, the residual sum of squares ($\ensuremath{\mathsf{RSS}\hspace*{0.2mm}}$) of a least-square estimation (we will impose conditions such that the estimator is well-defined). Hence, write
Also, for any element of $\boldsymbol{Z}$ indexed by $j\in[d_Z]$, and $\delta\in[0,1]$ define
The optimum splitting point $\delta^*$ of the index set $\mathcal{I}$ along direction $j$ is given by
Note that we might take $\delta^*\in\{Z_{i,j}:i\in\mathcal{I}\}$. Finally, define the left and right child nodes of $\mathcal{I}$ by \[ \ensuremath{\mathcal{I}}^{-}:=\{i\in\ensuremath{\mathcal{I}}:Z_{ij}\leq \delta^*\},\qquad\text{and}\qquad \ensuremath{\mathcal{I}}^+:=\{i\in\ensuremath{\mathcal{I}}:Z_{ij}> \delta^*\}. \]
Starting for an initial index set $\mathcal{B}\subseteq [n]$ and repeating the steps above, we end up with a partition of $[0,1]^{d_Z}$ into disjoint rectangles with axis-aligned sides (leaves). Furthermore, since the initial index set (root node) and the slipt direction $j$ are randomly chosen independent of the data, we denote by $\omega$ this independent source of randomness in the algorithm.
For $\boldsymbol{z}\in [0,1]^{d_Z}$, let $R(\boldsymbol{z},\omega)$ denotes the unique (random) leaf containing $\boldsymbol{z}$, i.e., $R(\boldsymbol{z},\omega)$ is a random rectangle that depends on the observation indexed by $\mathcal{B}$ and an external source of randomness. For an nonempty $\ensuremath{\mathcal{A}}\subseteq [n]$ such that $\ensuremath{\mathcal{A}}\cap\ensuremath{\mathcal{B}}=\emptyset$, define $\ensuremath{\mathcal{A}}(\boldsymbol{z},\omega)=\{i\in\mathcal{A}:\boldsymbol{Z}_i\in R(\boldsymbol{z},\omega)\}$, $\widehat{m}(\boldsymbol{x},\boldsymbol{z},\omega):=\boldsymbol{x}^{\ensuremath{\mathsf{T}}}\widehat{\boldsymbol\beta}(\boldsymbol{z},\omega)$, and $\widehat{\boldsymbol\beta}(\boldsymbol{z},\omega):= \widehat{\boldsymbol\beta}(\mathcal{A}(\boldsymbol{z},\omega))$.
More explicitly, write
where $\widehat{\boldsymbol{\Omega}}(\boldsymbol{z},\omega):= \frac{1}{|\ensuremath{\mathcal{A}}(\boldsymbol{z},\omega)|}\sum_{i\in\ensuremath{\mathcal{A}}(\boldsymbol{z},\omega)}\boldsymbol{X}_i\boldsymbol{X}_i^{\ensuremath{\mathsf{T}}}$ and $\widehat{\boldsymbol{\gamma}}(\boldsymbol{z},\omega):=\frac{1}{|\ensuremath{\mathcal{A}}(\boldsymbol{z},\omega)|}\sum_{i\in\ensuremath{\mathcal{A}}(\boldsymbol{z},\omega)}\boldsymbol{X}_iY_i$. Finally, for $B\geq 1$, we propose to estimate $m(\boldsymbol{x},\boldsymbol{z})$ using $\overline{m}(\boldsymbol{x},\boldsymbol{z}):= \frac{1}{B}\sum_{b=1}^B\widehat{m}(\boldsymbol{X},\boldsymbol{z},\omega_b) = \boldsymbol{X}^{\ensuremath{\mathsf{T}}}\overline{\boldsymbol\beta}(\boldsymbol{z})$ with
where $\{\omega_b:b\in[B]\}$ is an independent and identically sequence independent of $\{(Y_i,\boldsymbol{X}_i,\boldsymbol{Z}_i):i\in[n]\}$. We also investigate the properties of the following estimator: \[ \widecheck{\boldsymbol\beta}(\boldsymbol{z}):=\overline{\boldsymbol{\Omega}}(\boldsymbol{z})^{-1}\overline{\boldsymbol{\gamma}}(\boldsymbol{z}), \] where $\overline{\boldsymbol{\Omega}}(\boldsymbol{z}):= \frac{1}{B}\sum_{b=1}^B \widehat{\boldsymbol{\Omega}}(\boldsymbol{z},\omega_b)$ and $\overline{\boldsymbol{\gamma}}(\boldsymbol{z}):= \frac{1}{B}\sum_{b=1}^B \widehat{\boldsymbol{\gamma}}(\boldsymbol{z},\omega_b)$.
In this section, we state the main results of the paper and their underlying assumptions.
In Assumption (ref)(e), we require the minimum of observations per leaf to increase with the same size to ensure that each tree is consistent, as opposed to the consistency of the random forest. Specifically, we need the variance of $\widehat{\boldsymbol{\Omega}}(\boldsymbol{z})$ to vanish on each tree as the sample size increases. At the same time, $k$ must grow slower than the subsampling rate $s$ so that each cell is split enough times to make its diameter shrink toward zero, and the tree bias vanishes as the sample size increases.
Recall that, for the tree constructed with randomness $\omega$, $R(\boldsymbol{z},\omega)$ denotes the leaf containing $\boldsymbol{z}$. Define the map $\widetilde{\boldsymbol{\beta}}:[0,1]^d\to \R^{d_X}$ by \[ \widetilde{\boldsymbol\beta}(\boldsymbol{z}) := \nabla_x\E[Y|\boldsymbol{X}=\boldsymbol{x},\boldsymbol{Z}\in R(\boldsymbol{z},\omega)] = \E[\boldsymbol\beta(\boldsymbol{Z})|\boldsymbol{Z}\in R(\boldsymbol{z},\omega)];\qquad \boldsymbol{z}\in[0,1]^d. \] Note that in general we expect $\widetilde{\boldsymbol{\beta}}(\cdot)\neq \boldsymbol{\beta}(\cdot)$. However, due to the sample split we have that $\widehat{\boldsymbol\beta}(\cdot,\omega)$ is an unbiased estimator for $\widetilde{\boldsymbol\beta}(\cdot)$ even when $\widehat{\boldsymbol\beta}(\cdot,\omega)$ is biased for $\boldsymbol{\beta}(\cdot)$ (and $\widetilde{\boldsymbol\beta}(\cdot)$).
Although a tree construct using Algorithm (ref) is indeed honest and symmetric, there is no guarantee to be a $k$-PNN predictor. To see that, fix a test point $\boldsymbol{z}\in[0,1]^d$ and let $R(\boldsymbol{z},\omega)$ denote the unique leaf containing $z$. Recall that $R(\boldsymbol{z},\omega)$ is independent of $\mathcal{A}$-sample due to honesty. Then the number of observations in the $R(\boldsymbol{z},\omega)$ wich we denoted by $|\mathcal{A}(\omega,z)|$ follows a Binomial distribution conditional on $R(\boldsymbol{z},\omega)$ with $s$ trial an probability of success $p(\boldsymbol{z},\omega) :=\ensuremath{\mathbb{P}}\big(\boldsymbol{Z}\in R(\boldsymbol{z},\omega)|R(\boldsymbol{z},\omega)\big)$.
Therefore, the expected number of $\mathcal{A}$-sample observations in $R(\boldsymbol{z},\omega)$ conditional on $R(\boldsymbol{z},\omega)$ is given by $s p(\boldsymbol{z},\omega)$ while the number of $\mathcal{B}$-sample in $R(\boldsymbol{z},\omega)$ is between $k$ and $2k-1$ by construction. The question becomes how $|\mathcal{A}(\boldsymbol{z},\omega)|$ relates to $ k$. As it is shown in Lemma (ref), $|\mathcal{A}(\omega,z)|$ can be upper and lower bound in probability as \[ \ensuremath{\mathbb{P}}\left[s\left(\frac{s}{k}\right)^{-1/K(\alpha)}\lesssim|\mathcal{A}(\boldsymbol{z},\omega)|\lesssim s\left(\frac{s}{k}\right)^{-K(\alpha)}\right]\gtrsim 1 - \frac{1}{s(s/k)^{1/K(\alpha)}}. \] Set $k\asymp s^\eta$ for $\eta\in [0,1)$. For $\eta\in (1-K(\alpha),1)$ we have that $|A(\boldsymbol{z},\omega)|$ diverges as $s\to\infty$ with high probability. Precisely. \[ \ensuremath{\mathbb{P}}\left[s^{1-\frac{1-\eta}{K(\alpha)}}\lesssim|\mathcal{A}(\boldsymbol{z},\omega)|\lesssim s^{1-(1-\eta) K(\alpha)}\right]\to 1 \] For $\eta\in[0,1-K(\alpha)]$ and, in particular, $\eta=0$ (fixed $k$), the above bound is vacuous. It is not clear whether, for any fixed $k>0$, it is possible to claim that $|\mathcal{A}(\boldsymbol{z},\omega)|\gtrsim k$ with high probability.
We propose to estimate the covariance matrix $\boldsymbol\Sigma(\boldsymbol{z})$ appearing in Theorem (ref) using
where $\widehat{\boldsymbol{\Omega}}(\boldsymbol{z})$ is given by $\eqref{eq:beta_hat_definition}$. As for an estimator for $\widehat{\boldsymbol\Lambda}(\boldsymbol{z})$, define the residual function of the RF by $\boldsymbol{z}\mapsto \widehat{\epsilon}_i(\boldsymbol{z}):= Y_i - \boldsymbol{X}_i^{\ensuremath{\mathsf{T}}}\overline{\boldsymbol\beta}(\boldsymbol{z})$ for $i\in[n]$ and $\boldsymbol{z}\in[0,1]^{d_Z}$. Note any pair $\boldsymbol{z},\boldsymbol{z}'\in[0,1]^{d_Z}$ we observe $B$ independent realizations of $S(\boldsymbol{z},\boldsymbol{z}',\omega_b)$. So $\theta(\boldsymbol{z},\boldsymbol{z}')$ can be unbiased estimated by $\widehat{\theta}(\boldsymbol{z},\boldsymbol{z}'):=\frac{1}{B}\sum_{b=1}^B S(\boldsymbol{z},\boldsymbol{z}',\omega_b)$. So we propose to estimate $\widehat{\boldsymbol\Lambda}(\boldsymbol{z})$ using
Even when an observation $\boldsymbol{Z}_i$ is not used to grow a tree $\omega_b$, we can compute $S(\boldsymbol{z},\boldsymbol{Z}_i,\omega)$ by checking whether $\boldsymbol{Z}_i$ lies in $R(\boldsymbol{z},\omega)$ and if so count how many observations are in $R(\boldsymbol{z},\omega)$. Since we have $\E[\epsilon\theta(\boldsymbol{z},\boldsymbol{Z})X]=0$, we might consider the centered version of the estimator below given by
Suppose we wish to test that the conditional expectation does not depend on $\boldsymbol{Z}$ (homogeneous partial effect). Specifically, we are interested in the parametric null. \[ \mathcal{H}_0:\boldsymbol\beta(\boldsymbol{z}) = \boldsymbol\beta_0\qquad\text{for some $\boldsymbol\beta_0\in \R^{d_X}$ and all $\boldsymbol{z}\in[0,1]^d$} \] against the non-parametric alternative hypothesis $\mathcal{H}_1:\boldsymbol\beta(\boldsymbol{z})$ is Lipschitz. Following jFcZjZ2001; we propose to test $\mathcal{H}_0$ using the Generalized LRT, which is given as (after taking logs) \[ \Lambda(\mathcal{H}_0):=\frac{n}{2}\log\frac{\ensuremath{\mathsf{RSS}\hspace*{0.2mm}}_0}{\ensuremath{\mathsf{RSS}\hspace*{0.2mm}}}, \] where $\ensuremath{\mathsf{RSS}\hspace*{0.2mm}}_0 := \sum_{i=1}^n (Y_i - \boldsymbol{X}_i^{\ensuremath{\mathsf{T}}}\widetilde{\boldsymbol\beta})^2 $, $\ensuremath{\mathsf{RSS}\hspace*{0.2mm}} := \sum_{i=1}^n \left[Y_i - \boldsymbol{X}_i^{\ensuremath{\mathsf{T}}}\overline{\boldsymbol\beta}(\boldsymbol{Z}_i)\right]^2 $ and $\widetilde{\boldsymbol\beta}$ is the OLS estimator.
The idea is to explore the fact that for a Lagrange Multiplier (LM) type test, we are only required to estimate the model under the null and obtain a consistent estimator under both the null and alternative. When the null completely characterizes $\boldsymbol\beta(\boldsymbol{z})$ or when $\boldsymbol\beta(\boldsymbol{z})$ is known up to finite unknowns (parametric), we propose a somewhat canonical test.
Under $\mathcal{H}_0:\boldsymbol\beta(\boldsymbol{z})$ is not a function of $z$, we have that $\E[\boldsymbol{Z}(Y-\boldsymbol{X}^{\ensuremath{\mathsf{T}}}\widetilde{\boldsymbol\beta})=0$ for some unknown $\widetilde{\boldsymbol\beta}\in\R^{d_X}$. Let $\widehat{\epsilon}_{i}^ {OLS}:=Y_i - \boldsymbol{X}_i^{\ensuremath{\mathsf{T}}}\widehat{\boldsymbol\beta}_{OLS}$ for $i\in[b]$ where $\widehat{\boldsymbol\beta}_{OLS}$ is the OLS estimator of $Y$ regressed onto $X$
So, we can use the sample moment condition below as a basis for constructing our test statistic. \[ M=\sum_{i=1}^n \widehat{\epsilon}_{i}^ {OLS}\boldsymbol{Z}_i, \] because under $\mathcal{H}_0$ \[ M/\sqrt{n} = L\frac{1}{\sqrt{n}}\sum_{i=1}^n \epsilon_iW_i + o_\ensuremath{\mathbb{P}}(1);\qquad W_i := (\boldsymbol{Z}_i^{\ensuremath{\mathsf{T}}}, \boldsymbol{X}_i^ {\ensuremath{\mathsf{T}}})^{\ensuremath{\mathsf{T}}} \] where $L:=(I_d: -\E[Z\boldsymbol{X}^{T}]\big(\E[\boldsymbol{X}\boldsymbol{X}^{T}\big)^{-1})$ is an $(d_Z (d_Z+d_X))$ matrix. Since $\frac{1}{\sqrt{n}}\sum_{i=1}^n \epsilon_iW_i\overset{d}{\longrightarrow} \mathsf{N}(0,\E[\epsilon^2WW^ {\ensuremath{\mathsf{T}}}])$, we have that $M/\sqrt{n}\overset{d}{\longrightarrow} \mathsf{N}(0,L\E[\epsilon^2WW^ {\ensuremath{\mathsf{T}}}]L^{\ensuremath{\mathsf{T}}})$ therefore \[ M^ {\ensuremath{\mathsf{T}}} (nV)^{-1} M\overset{d}{\longrightarrow}\chi^2_{d_Z};\qquad V:=L\E[\epsilon^2WW^ {\ensuremath{\mathsf{T}}}]L^{\ensuremath{\mathsf{T}}} \] Let $\widehat{\epsilon}_{i}^{RF}:=Y_i - \boldsymbol{X}_i^{\ensuremath{\mathsf{T}}}\overline{\boldsymbol\beta}$ be the residuals of the RF. Under the alternative (and the null) $\overline{\boldsymbol\beta}(\boldsymbol{z})\overset{\ensuremath{\mathbb{P}}}{\longrightarrow} \boldsymbol\beta(\boldsymbol{z})$ then $\epsilon_{i}^{RF}\overset{\ensuremath{\mathbb{P}}}{\longrightarrow}\epsilon_i$ for $i\in[n]$. Define the plug-in estimators \[ \widehat{L}:=\left[\boldsymbol{I}_d: -\sum_{i=1}^nZ_i\boldsymbol{X}^{\ensuremath{\mathsf{T}}}_i\left(\sum_{i=1}^n\boldsymbol{X}_i\boldsymbol{X}^{\ensuremath{\mathsf{T}}}_i\right)^{-1}\right],\qquad \widehat{V}:=\widehat{L}\frac{1}{n}\sum_{i=1}^ n\left[(\widehat{\epsilon}_{i}^{RF})^2W_iW^{\ensuremath{\mathsf{T}}}_i\right]\widehat{L}^{\ensuremath{\mathsf{T}}} \]
Hence, the test statistics become \[ T:= M^ {\ensuremath{\mathsf{T}}} (n \widehat{V})^{-1} M\overset{d}{\longrightarrow}\chi^2_{d_Z} \] as $n\to\infty$ under $\mathcal{H}_0$.
Note that the same process works to test $\mathcal{H}_0:\boldsymbol\beta(\boldsymbol{z})=\boldsymbol\beta_0(\boldsymbol{z})$ for some known function $\boldsymbol{z}\mapsto\boldsymbol\beta_0(\boldsymbol{z})$.
Partition the control variables as $\boldsymbol{Z}= (\boldsymbol{Z}',\boldsymbol{Z}'')$ where $\boldsymbol{Z}''= (Z_{1}'',\dots Z_{d_{\boldsymbol{Z}''}}'')$ and $\boldsymbol{Z}''_j$ are discrete random variables with $m_j\geq 2$ categories for $j\in[d_{\boldsymbol{Z}''}]$. Without loss of generality we may assume that $Z_{j}''$ is supported on $\mathcal{S}_j=\{0,1/(m_j-1),1/(m_j-2), \dots, 1\}$ and thus $\boldsymbol{Z}''$ has support on $\mathcal{S} = \bigtimes_{j=1}^{d_{\boldsymbol{Z}''}} \mathcal{S}_j$.
Let $R=\bigtimes_{j=1}^{d_Z}[a_j,b_j]$ denote a rectangle in $[0,1]^{d_Z}$. For convenience, we can define the length of a rectangle along a discrete variable $\boldsymbol{Z}''_j$ as the fraction of the categories in the rectangle. Specifically $0\leq \operatorname{diam}(R)_j := (|\mathcal{S}_j\cap [a_j,b_j]|-1)/(m_j-1)\leq 1$.
When a discrete variable $Z_j$ is chosen to split we replace $\eqref{eq:MSE_continuos}$ by the condition
and identify the left and right child nodes by the index sets
Fix a test point $\boldsymbol{z}=(\boldsymbol{z}',\boldsymbol{z}'' )\in[0,1]^{d_{\boldsymbol{Z}'}}\times \mathcal{S}$ and let $M_j(\boldsymbol{z})$ denote the number of splits along the $j$-th discrete variable to form the (unique) leaf containing $z$. Since $s/k\to \infty$ we have that some large $s$ \[ \ensuremath{\mathbb{P}}(M_j(\boldsymbol{x})< m_j-1)\leq \ensuremath{\mathbb{P}}\left( M_{j}(\boldsymbol{z})\leq \tfrac{(1-\delta)\pi}{d_Z} \frac{\log(s/(2k-1))}{\log(1/(1-\alpha))}\right)\leq \left(\frac{s}{2k-1}\right)^{-\frac{\delta^2\pi^2}{2d_Z^2\log(1/(1-\alpha))}}. \] Let $\mathcal{E}_D(\boldsymbol{x}) = \bigcap_{j\in[d_Z'']}\{M_j(\boldsymbol{x}) = m_j-1\}$, then by the union bound we have \[ \ensuremath{\mathbb{P}}(\mathcal{E}_D(\boldsymbol{x})) \gtrsim 1 -d_Z''\left(\frac{s}{k}\right)^{-\frac{\delta^2\pi^2}{2d_Z^2\log(1/(1-\alpha))}} \] i.e., all discrete variables have a single class in the leaf $R(\boldsymbol{z},\omega)$ with probability at least $1-\left(\frac{s}{k}\right)^\frac{-\delta^2\pi^2}{2d_Z^2\log(1/(1-\alpha))}$. Hence, the bias on the leaf $R(\boldsymbol{z},\omega)$ can be upper bounded on $\mathcal{E}_D(\boldsymbol{x})$ as \[ \max_{(u',u'')\in R((z',z''),\omega)}\|\boldsymbol\beta(u',u'')-\boldsymbol\beta(z',z'')\|\leq \|\boldsymbol\beta(u',z'')-\boldsymbol\beta(z',z'')\|\leq C\operatorname{diam}(R(\boldsymbol{z},\omega)) \]
We conducted a simulation study to test the validity of our results in finite sample. We consider several specifications with different sample sizes. Specifically, we set the number of observations to be $n\in\{250,500,1000\}$. The dimension $d_Z$ of the variables determining model heterogeneity is $d_Z=\{1,2,3,5\}$. Each simulated model is estimated by the Random Forest as per the Algorithm (ref) with $B=3000$ subsampling replications. We consider a fraction of $s=0.8$ observations in each subsample, $k=s^{1/6}$ and $\alpha=0.005$. The number of Monte Carlo replications is 500.
The simulated models are determined by the general equation
where $U_i$ is an independent and normally distributed random variable with zero mean and standard deviation equal to 0.5. The scalar continuous treatment variable is $X_i\sim\textsf{Uniform}(0,1)$. $\boldsymbol{Z}_i$ is a random vector of $d_Z$ mutually independent random variables uniformly distributed between 0 and 1.
We set $\beta_0(Z_i)=0$ for all $i=1,\ldots,n$, and consider the following specifications for the slope parameter $\beta_1(\boldsymbol{Z}_i)$:
The results are reported in Figure (ref) and Tables (ref)--(ref). Figure (ref) reports results for Model I where $p=1$. It illustrates the median slope estimation across the Monte Carlo simulations, as well as the 95% confidence bands. It is evident the estimation is precise in this case.
Table (ref) presents the averages derived from the Monte Carlo simulations for various descriptive statistics pertinent to goodness-of-fit assessments. The evaluated statistics include the mean, standard deviation, kurtosis, and skewness. For Model I, Panel (I.a) examines the estimated residuals from the model fit, Panel (I.b) addresses the estimation error associated with the varying intercept of the model, and Panel (I.c) pertains to the estimation error for the varying slope coefficient. The subsequent panels replicate these findings for Models II, III, and IV.
Several facts emerge from the table. First, the model approximation is satisfactory even in samples as small as 250 observations. Notably, the average of the residual standard deviation is close to 0.5, which is the true value. Furthermore, the average estimated kurtosis is approximately three, and the average estimated skewness is close to zero, indicating that the residuals are approximately normally distributed. Additionally, the results suggest that the estimation of the varying intercept is more precise than the slope estimation, as the average standard deviation of the errors is significantly smaller for the former than for the latter. Finally, as expected, the performance of the method slightly deteriorates as the dimension of $\boldsymbol{Z}$ increases.
Tables (ref) and (ref) present results for several test values pertinent to the vector $\boldsymbol{Z}$. We analyze 11 points where $\boldsymbol{Z}$ constitutes a $5 1$ that spans from a vector of zeros to a vector of ones, with $0.1$ increments. For simplicity and without loss of generality, we assume that all elements of the test points are equal. Table (ref) reports the average bias and the mean squared error (MSE) for the varying intercept coefficient evaluated at each test point for different sample sizes. It is evident that the intercept estimation is precise and improves with increasing sample size. Table (ref) reports the average bias and the mean squared error (MSE) for the varying slope coefficient evaluated at each test point for different sample sizes. The results in the table corroborate our previous conclusion that as the sample size increases, the method's performance improves. Finally, there is evidence that the function approximation is better for points closer to the center of the distribution of $\boldsymbol{Z}$.
We present coverage results in Tables (ref) and (ref). Specifically, we provide the average 90% and 95% coverage across Monte Carlo simulations for both the intercept and the slope parameters. Table (ref) details the coverage for the intercept, while (ref) presents the corresponding results for the slope parameter. It is evident from the tables that the coverage probabilities for the intercept are close to the expected values, even with small sample sizes and an increased number of variables. This finding supports our previous simulation results, indicating minimal bias in intercept estimation. Conversely, the coverage for the slope parameter is significantly underestimated, particularly for values that lie farther from the center of the covariate distribution. This observation aligns with the larger biases reported in Table (ref). It is important to mention that such narrow coverage probabilities have also been reported in the simulations in rFjTsAsW2021.
Finally, to evaluate the finite sample performance of the Lagrange Multiplier homogeneity (linearity) test described in Section (ref). We simulated linear (homogeneous) models with $p\in\{1,2,3,5\}$ covariates. The results are reported in Figure (ref), which illustrates the size discrepancy relative to the nominal size. As observed, the size distortions remain negligible.
In this section, we illustrate our proposed methodology utilizing two distinct datasets. The first is a synthetic dataset generated by a widely used model specifically designed to evaluate nonparametric methods. The second is a real dataset concerning the economic convergence of Brazilian municipalities.
The first dataset consists of observations generated according to the following model:
where $\boldsymbol{Z}_i$ is a vector of mutually independent uniform random variables taking values on $[0,1]^5$, and $U_i$ is a zero-mean normally distributed random variable with unit variance. The model for $\beta(\boldsymbol{Z}_i)$ has been widely employed in the literature to evaluate semi-parametric models. In this context, heterogeneity is jointly influenced by interactions, quadratic forms, and a robust linear signal. Figure (ref) illustrates the three sources of heterogeneity. See, for example, rFjTsAsW2021 for a similar data-generating process.
To evaluate the performance of the model presented in this paper, we consider various sample sizes (2000, 4000, 8000, 16000, 32000). Figure (ref) illustrates the empirical distribution of $\boldsymbol\beta(\boldsymbol Z_i)$, $i=1, \ldots, 32000$. We estimate the locally linear random forest model using the same hyperparameter settings as those utilized in the Monte Carlo simulation.
Figure (ref) shows results concerning the estimation parameters. Panels (a) and (b) illustrate the scatter plot of the fitted $\boldsymbol\beta(\boldsymbol{Z})$ against the true values for $n=2000$ and $n=36000$, respectively. A linear regression line is also included. Figure (ref) reports the evolution of the bias and the MSE as a function of the sample size. As we can observe from both figures, the estimation improves as the sample size increases, which aligns with established statistical theory. However, as anticipated by our theoretical results, the convergence to the true values is notably slow.
We illustrate our methodology by testing heterogeneity in the convergence among Brazilian municipalities between 1970 and 2000. Our starting point is the simplified convergence equation presented in rjBxS1992:
where $Y_{i,t}$ is the per capita income of region $i$ in period $t$, $a_{i}$ is associated with the steady-state level of per capita income and the rate of technological progress, $\phi_{i}$ is a parameter related to the time trend determined by the technological progress, and $U_{i,t}$ is the random term. Convergence corresponds to the parameter $\gamma_i$.
From a conceptual point of view, two alternative assumptions determine the most important distinction of convergence concepts. First, we can assume that $a_{i}=a$, $\gamma_i=\gamma$, and $\phi_{i}=\phi $, i.e., that the basic parameters of preference and technology are the same for all economies represented in the sample. This is when $\gamma<0$ represents unconditional convergence - a situation where poorer regions tend unconditionally to grow more quickly than richer ones. Alternatively, we can state a weaker assumption, allowing for possible differences in the steady state across the economies considered and heterogeneity in the convergence rate. In terms of equation ((ref)), $a_{i}$, $\gamma_i$ and $\phi _{i}$ are allowed to vary across different regions.
We estimate ((ref)) in a cross-section setup, where there is no identifiable time trend, and we are not able to distinguish between $\phi_{i}$ and $a_{i}$. Thus, we estimate the following equation:
The data originate from the Brazilian Demographic Censuses conducted in 1970 and 2000. The geographical units were adjusted to account for the reorganization of Brazilian territory throughout this time frame. In 1970, Brazil was composed of 3,951 municipalities. By 2000, the count had increased to 5,507 municipalities. Consequently, all data collected in 2000 were aggregated to align with the municipal structure as it existed in 1970. Our dependent variable is defined as the average growth in per capita income from 1970 to 2000 for each municipality. The independent variable utilized in this analysis is the logarithm of the per capita income level recorded in 1970.
Figure (ref) illustrates the differences across municipalities in terms of growth rate. Figure (ref) shows that the variations in growth do not coincide with the administrative state frontiers. There are substantial variations within many of the Brazilian states.
We employ the semi-parametric approach outlined in the previous section to assess conditional convergence. Variations in preferences and technological parameters are endogenously incorporated into the analysis through geographical proximities. The underlying assumption suggests that cities situated in close proximity experience similar steady states. Within our modeling framework, we estimate ((ref)), recognizing that $\alpha _{i}$ represents a semi-parametric function of the latitude and longitude coordinates. The results are displayed in Figures (ref) and (ref), which report the heterogeneity patterns in the intercept and the convergence parameter. The estimated geographic heterogeneity varies remarkably, with notable clusters of high-income municipalities in the central and southern parts of the country. Significant geographic differences across municipalities do not coincide with the geographic structure of the Brazilian states.
This paper presents a semi-parametric framework for estimating heterogeneous partial effects, combining the flexibility of machine learning techniques with the interpretability of traditional parametric models. By utilizing a Random Forest-based methodology, we provide a robust and adaptable approach to capturing complex heterogeneity in the relationship between explanatory variables and outcomes.
Our theoretical contributions establish key consistency and asymptotic normality results, ensuring the reliability of our estimator. Importantly, our framework accommodates both continuous and discrete covariates while maintaining desirable statistical properties. The Monte Carlo simulations demonstrate the method’s accuracy, even in moderate sample sizes, and highlight the precision of our approach in recovering varying intercepts and slopes. Additionally, our empirical analysis of Brazilian municipal economic convergence underscores the practical relevance of our method, revealing substantial geographic heterogeneity in growth dynamics.
Future research can extend this methodology to settings with dependent data and high-dimensional covariates, broadening its applicability to dynamic panel models and network data. The current approach can accommodate high-dimensional controls ($\boldsymbol{Z}$) under additional conditions on the underlying data-generating process. Specifically, it would require that (i) only a small number of controls are relevant to the model (sparsity assumption) and (ii) the relevant controls are independent of the irrelevant controls and the outcome variable. The latter is a strong assumption in most applications, so we did not pursue this route here.
Regarding the dependent data across observations, the subsampling and sample split steps must be adjusted to preserve both the data dependence structure and the independence of the two samples. A block-subsampling technique is a natural candidate; however, implementing this would significantly alter the proof techniques for the subsequent steps and might obscure the main idea of the paper. Furthermore, refinements in inference procedures for heterogeneous effects, including the construction of a uniformly valid confidence interval, would enhance the model’s applicability in applied research. Additionally, considering Lemma (ref)(d), a significant modification to the algorithm would be necessary to achieve uniform convergence of the proposed estimators.
Overall, this paper offers a flexible, interpretable, and computationally efficient tool for studying heterogeneity in the partial effect of a variable of interest, bridging the gap between parametric and nonparametric estimation.