Bayesian Hierarchical Random Forests

• Bayesian Hierarchical Random Forests are ensemble models that merge the flexibility of decision trees with a Bayesian treatment of uncertainty and structured regularization.
• They use a two-level approach where stable top-level splits are fixed while deeper branches randomize, enhancing predictive accuracy and reducing computational costs.
• Applications include high-dimensional regression and spatial forecasting, where hierarchical shrinkage and integrated CAR models improve inference and model performance.

Bayesian Hierarchical Random Forests (BHRF) denote a class of statistical machine learning models that merge the flexibility of decision-tree ensembles—principally Random Forests—with a principled Bayesian hierarchical treatment of uncertainty, regularization, and, when necessary, domain-specific covariance structure. These approaches are motivated by the need to combine interpretable, nonlinear predictive power with robust uncertainty quantification, shrinkage for high-dimensional feature spaces, and the capacity to model complex dependencies such as spatial autocorrelation or hierarchical group structures.

1. Bayesian Nonparametric Ensemble Frameworks

Bayesian forests extend random forests by treating tree ensembles as draws from a nonparametric Bayesian model. In such frameworks, the data-generating process (DGP) is modeled as discrete with an unknown probability law. The Bayesian bootstrap and Dirichlet priors induce a posterior over finite-support distributions, and subsequently, a distribution over tree structures via the mechanism of weighted CART splits. Specifically, let $D=\{(x_i,y_i)\}_{i=1}^n$ be the observed data. The nonparametric prior places weights $\omega_i$ (drawn from a limiting Dirichlet process) on each observation, which are implemented as $\theta_i\sim \mathrm{Exp}(1)$, $\omega_i = \theta_i / \sum_{j=1}^n \theta_j$. Each tree $\mathcal{T}(\theta)$ in the Bayesian Forest is the CART tree fitted on this weighted sample. The ensemble posterior arises from repeated independent resampling, yielding both predictions and a quantification of model uncertainty (Taddy et al., 2015).

2. Hierarchical Character of Bayesian Forests

BHRF models exhibit a natural hierarchy in the inferred structure of ensembles. Statistical analysis shows top-level splits are extremely stable as node size increases: for a node with $n$ observations, the probability that the best split (by impurity minimization) differs from the sample CART version decays as $p/\sqrt{n}\, e^{-n}$, where $p$ is the number of features. High-level (shallow) nodes are "clamped" and show little to no posterior uncertainty, while lower branches exhibit greater variability. This property motivates empirical Bayesian strategies: fixing the trunk at its modal value and only randomizing the lower-level branches ("Empirical Bayesian Forests," EBFs). This modification maintains predictive accuracy while substantially reducing computational cost (Taddy et al., 2015).

3. Hierarchical Shrinkage and Regularization of Tree Rules

Advanced BHRF formulations exploit hierarchical priors to perform aggressive regularization, especially when ensembles generate large numbers of indicator rules, many of which correspond to noise. This is operationalized in "Tree Ensembles with Rule Structured Horseshoe Regularization" (Nalenz et al., 2017), where the model is defined as: $$y_i \mid \beta, \alpha, \sigma^2 \sim \mathcal{N}(\mu_i, \sigma^2), \quad \mu_i = \sum_{j=1}^J R_j(x_i)\beta_j + \sum_{k=1}^K x_{i,k}\alpha_k$$ with $R_j(x)$ indicating activation of rule $j$ (extracted from both random forest and gradient-boosted trees). Regression weights are assigned a three-layer hierarchical horseshoe prior: $$\theta_m \mid \lambda_m, \tau, \sigma^2 \sim \mathcal{N}(0, \lambda_m^2 \tau^2 \sigma^2)$$

$\lambda_m \sim \mathcal{C}^+(0, A_m)$

$\tau \sim \mathcal{C}^+(0,1)$

with $A_m$ encoding rule-specific shrinkage based on rule complexity and empirical prevalence. Noise rules are shrunk close to zero, controlling overfitting even with thousands of candidate rules (Nalenz et al., 2017).

4. Hierarchical Bayesian Spatial Random Forests

Bayesian hierarchical random forests are adapted for spatial and areal prediction problems via hybridization with Bayesian conditional autoregressive (CAR) spatial models. The CAR-Forest architecture (MacBride et al., 2023) fuses nonparametric random forest predictors with Bayesian spatial random effects. The algorithm proceeds by alternating between (A) random forest fitting, which captures complex nonlinear mean structure, and (B) Bayesian CAR fitting on the residuals, which accounts for spatial autocorrelation: $$Y_k \mid \beta_0, \phi_k, \sigma^2, \widehat{m}^{(r)}_k \sim \mathcal{N}(\beta_0 + \widehat{m}^{(r)}_k + \phi_k, \sigma^2)$$

$$\phi \sim \mathcal{N}\left(0, \tau^{-1} Q(\rho)^{-1}\right)$$

where $Q(\rho)$ encodes spatial precision via neighborhood structure and a smoothing parameter $\rho$. Iterative updates allow the model to partition variation between covariate-driven and latent spatial components (MacBride et al., 2023).

5. Posterior Inference, Implementation, and Interpretation

Posterior inference in BHRF typically relies on repeated draws, either via direct resampling (in the Bayesian Forest/EBF setting) or via efficient Gibbs samplers (for models with explicit hierarchical structure). In the "horseshoe regularization" context, all update steps admit conjugate forms, enabling scalable MCMC. Posterior summaries include average predictions, credible intervals, and variable or rule importance metrics computed from the posterior draws: $$I(R_j) = |\beta_j|\,\mathrm{sd}(R_j(x)), \quad I(x_k) = |\alpha_k|\,\mathrm{sd}(x_k)$$

$$J(x_k) = I(x_k) + \sum_{j: k\in Q_j} \frac{1}{|Q_j|}\,I(R_j)$$

Graphical summaries such as RuleHeat plots and marginal effect plots, as well as Shapley-value analyses for interpretability, are implemented in dedicated R packages (e.g., "horserule") (Nalenz et al., 2017).

6. Empirical Results and Applications

Empirical evaluation demonstrates that BHRF models can achieve state-of-the-art performance in various settings. On high-dimensional regression and classification tasks, aggressive regularization via hierarchical horseshoe priors yields improvements over RuleFit, BART, and random forests (Nalenz et al., 2017). In spatial prediction settings, the CAR-Forest model outperforms CAR, random forest, and geographically weighted random forest, delivering reduced RMSE, improved coverage, and narrower predictive intervals for small-area estimation (e.g., Scottish housing prices: RMSE = £40,303, MAE = £16,623, 95% coverage = 0.944, and interval width = £140,875) (MacBride et al., 2023). EBFs retain predictive performance close to full Bayesian forests while offering large computational savings, especially in massive-data regimes (Taddy et al., 2015).

7. Extensions and Theoretical Considerations

BHRF frameworks are extensible to multivariate and spatio-temporal domains, accommodating separate forests for each margin and joint CAR models for complex dependencies. Split allocation between forest and explicit model components can be exploited for interpretability and confounder-control in ecological regression applications. The theoretical analysis reveals that as node size increases, the top-level splits in a Bayesian decision tree ensemble concentrate superexponentially, justifying empirical Bayesian approximations and modal trunk fixing. When the true relationship is linear, supplementing forests with explicit linear fixed effects in the hierarchical model is recommended to maintain efficiency (Taddy et al., 2015, Nalenz et al., 2017, MacBride et al., 2023).