Wednesday, October 7, 2026

Introduction to the Random Forest algorithm (for Hydrologists)

Random forests are among the most used machine-learning tools in hydrology, and among the least understood by the people who use them. This post is a guided reading of how they work, written for hydrologists who want the logic before the library call.

I follow one source closely: Chapter 8, Tree-Based Methods, of An Introduction to Statistical Learning (ISL) by Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani, together with the Stanford video lectures in which Hastie and Tibshirani teach that chapter. The book exists in an R edition (2nd ed., 2021) and a Python edition (ISLP, 2023, with Jonathan Taylor); both are free as PDFs from statlearning.com. The chapter numbering and notation below are theirs. Where I depart from them, it is to add a hydrological reading, and I say so. The relevant lectures are embedded section by section; the whole Chapter 8 series is in this YouTube playlist.

The chapter builds in four steps, and so does this post: a single tree, then many trees averaged (bagging), then many decorrelated trees averaged (random forests), then trees grown sequentially (boosting). Each step fixes a weakness of the previous one.


1. A single regression tree (ISL §8.1)

A regression tree cuts the predictor space into \(J\) non-overlapping boxes \(R_1, \dots, R_J\) and predicts, in each box, the mean of the training responses that fall in it. Think of predicting mean annual runoff from basin area, mean elevation, aridity index and forest cover: the tree asks is aridity < 0.8?, then is elevation > 1500 m?, and so on, until each leaf holds a group of similar basins.

The boxes are chosen to minimise the residual sum of squares:

\[ \mathrm{RSS} = \sum_{j=1}^{J} \sum_{i \in R_j} \left( y_i - \hat{y}_{R_j} \right)^2 \]

Finding the best partition is computationally infeasible, so the tree is grown by recursive binary splitting: at each step, take the single predictor \(X_j\) and cutpoint \(s\) that most reduce RSS, split, and repeat inside each half. It is greedy and top-down; it never looks back to see whether an earlier split was a poor choice.

Hastie and Tibshirani, Decision Trees (Stanford Statistical Learning, Ch. 8, 14:37).

A tree grown until each leaf holds a handful of points fits the training data almost perfectly and predicts new data badly. ISL's remedy is cost-complexity pruning: grow a large tree \(T_0\), then for each value of a penalty \(\alpha\) find the subtree \(T\) that minimises

\[ \sum_{m=1}^{|T|} \sum_{i:\, x_i \in R_m} \left( y_i - \hat{y}_{R_m} \right)^2 + \alpha |T| \]

where \(|T|\) is the number of leaves. \(\alpha\) plays the same role as \(\lambda\) in the lasso, and it is chosen by cross-validation.

Pruning a Decision Tree (Ch. 8, 11:45).

For classification the same machinery applies, with the Gini index or cross-entropy replacing RSS; the lecture Classification Trees and Comparison with Linear Models covers it.

Hastie and Tibshirani are frank about the balance sheet. Trees are easy to explain, mirror human decision-making, handle qualitative predictors without dummy variables, and can be drawn. But a single tree is usually less accurate than other methods, and it is non-robust: change the data a little and the tree can change a lot. In statistical language, trees have high variance. Everything that follows in the chapter is a way of killing that variance.

2. Bagging (ISL §8.2.1)

The idea rests on an elementary fact. If \(Z_1, \dots, Z_n\) are independent, each with variance \(\sigma^2\), their mean has variance \(\sigma^2/n\). Averaging reduces variance. If we had \(B\) independent training sets, we could grow \(B\) trees and average their predictions.

We do not have \(B\) training sets, so we fake them with the bootstrap (ISL Chapter 5): draw \(B\) samples of size \(n\), with replacement, from the one training set we have; grow a deep, unpruned tree on each; average:

\[ \hat{f}_{\mathrm{bag}}(x) = \frac{1}{B} \sum_{b=1}^{B} \hat{f}^{*b}(x) \]

For classification, average becomes majority vote. Each tree is deliberately overfitted (low bias, high variance); the averaging removes the variance. \(B\) is not a tuning parameter in the usual sense: a large \(B\) does not overfit, it simply costs time. You stop when the error has stabilised, often a few hundred trees.

Bagging also gives an error estimate for free. Each bootstrap sample leaves out, on average, about one third of the observations (the fraction tends to \(1/e \approx 0.368\)). For each observation, predict it using only the roughly \(B/3\) trees that did not see it. The resulting out-of-bag (OOB) error is a valid test-error estimate, and for large \(B\) it is essentially equivalent to leave-one-out cross-validation, as the authors note in the lecture below.

A caution the book does not need to make, but hydrologists do: the bootstrap assumes the observations are exchangeable. Daily streamflow is not. If neighbouring days land in both the bootstrap sample and the OOB set, the OOB error is optimistic. With time series, a block or year-based hold-out is the honest test.

3. Random forests (ISL §8.2.2)

Random forests improve bagging with one small change: at each split, the tree may only choose among a random subset of \(m\) of the \(p\) predictors, drawn afresh at every split. Typical choices are \(m \approx \sqrt{p}\) for classification and \(m \approx p/3\) for regression; \(m = p\) gives back bagging.

Bootstrap Aggregation (Bagging) and Random Forests (Ch. 8, 13:45): the core lecture for this post.

Why would forbidding a tree from seeing the best predictor help? Because bagged trees are correlated. Suppose one predictor is very strong (in a runoff problem, precipitation almost always is). Every bagged tree will then split on it first, and the trees will look alike. Averaging similar things does not reduce variance much. The arithmetic, given in The Elements of Statistical Learning (§15.2) and implicit in ISL, is the variance of an average of \(B\) identically distributed variables with pairwise correlation \(\rho\):

\[ \mathrm{Var}\!\left( \frac{1}{B}\sum_{b=1}^{B} \hat{f}^{*b} \right) = \rho\,\sigma^2 + \frac{1-\rho}{B}\,\sigma^2 \]

As \(B\) grows the second term vanishes, but the first does not. Bagging attacks the second term; random forests attack \(\rho\). By hiding the strong predictor from most splits, the other predictors get a chance, the trees diversify, and the average becomes more stable. In the lecture Hastie illustrates this on a high-dimensional gene-expression example, where \(m = \sqrt{p}\) clearly beats bagging; the gain is largest when many predictors are correlated.

The practical consequences:

  • Like bagging, a random forest does not overfit by adding trees; \(B\) only needs to be large enough.
  • \(m\) is the one parameter worth tuning, using the OOB error.
  • Trees are grown deep and unpruned. Minimum leaf size (typically 5 for regression) is a secondary control.
  • The price is interpretability: a forest of 500 trees cannot be drawn.

4. Variable importance, and how to misread it

To recover some interpretability, ISL records, for each predictor, the total decrease in RSS (or Gini index) produced by splits on it, averaged over all \(B\) trees. Ranked and plotted, this is the familiar variable-importance chart that ends half the hydrological ML papers I review. The lecture in §5 below discusses it alongside boosting.

Three warnings, the first from ISL's spirit and the others from the later literature:

  1. Importance is predictive, not causal. A high score says the forest found the variable useful for splitting this dataset. It does not say the variable drives the process.
  2. Correlated predictors share or steal credit. Elevation, temperature and snow fraction carry overlapping information; the forest spreads importance among them in ways that depend on \(m\) and on chance.
  3. Impurity-based importance is biased toward continuous predictors and categorical ones with many levels (Strobl et al., 2007). Permutation importance (Breiman, 2001), which measures how much the OOB error grows when one predictor's values are shuffled, is the safer default.

If the question is physical, the variable-importance plot is the beginning of an argument, not its conclusion.

5. The other branch: boosting and BART (ISL §8.2.3–8.2.4)

Chapter 8 does not stop at forests, and a reader should know why. Bagging and random forests grow trees independently and average them. Boosting grows them sequentially: each small tree is fitted to the residuals of the current model, and only a shrunken fraction of it is added:

\[ \hat{f}(x) \leftarrow \hat{f}(x) + \lambda\, \hat{f}^{b}(x), \qquad r_i \leftarrow r_i - \lambda\, \hat{f}^{b}(x_i) \]

The model learns slowly. It has three tuning parameters: the number of trees \(B\) (which, unlike in a forest, can overfit), the shrinkage \(\lambda\) (typically 0.01 or 0.001), and the tree depth \(d\) (often 1 or 2; with \(d = 1\), stumps, the model is additive). Boosting usually edges out random forests in accuracy but needs more care. XGBoost and LightGBM are its modern industrial descendants.

Boosting and Variable Importance (Ch. 8, 12:03).

The 2nd edition of the book adds Bayesian Additive Regression Trees (BART), which mixes both ideas: many trees, each repeatedly perturbed to fit the partial residual left by the others, with the output being a posterior distribution, hence uncertainty bands, at little tuning cost. BART is not in the 2014 video series.

6. A hydrologist's reading

What follows is mine, not the book's. Random forests are excellent at some hydrological jobs and structurally unable to do others, and the reasons follow directly from the construction above.

Where they shine. Regionalisation (predicting signatures or parameters in ungauged basins from catchment attributes); digital soil mapping; land-cover and snow-cover classification; gap-filling of station records; and post-processing, that is, learning the residuals of a process-based model such as one built in GEOframe. In all these cases the predictors are tabular, the relations are non-linear with interactions, and the forest needs almost no tuning.

Where they fail by design.

  • No extrapolation. A forest predicts averages of training responses, so its output is bounded by what it has seen. A flood larger than any in the record, or a climate warmer than the calibration period, is outside its reach. This is the most important single fact about trees for climate-impact work.
  • No conservation laws. Nothing ties predicted runoff, evapotranspiration and storage change to a closed water budget. Each output is fitted separately.
  • No memory. A forest maps inputs to output; dynamics must be put in by hand through lagged predictors, antecedent indices and the like.
  • Dependent data. As noted for OOB error, spatial and temporal autocorrelation make random cross-validation optimistic. Use spatially blocked or leave-basins-out schemes.

Uncertainty. A plain forest gives a point estimate. Quantile regression forests (Meinshausen, 2006) keep the full distribution of responses in each leaf and return predictive quantiles at no extra cost, which is often what a hydrological application actually needs.

My own conclusion is that forests are best used with physics, not instead of it: as a fast, honest benchmark that a process model must beat, and as a tool for learning what the process model gets wrong.

7. Try it: a synthetic Budyko test

The example below generates 500 synthetic basins whose runoff follows Fu's form of the Budyko curve, plus a small forest effect and noise. Elevation is a decoy with no effect. It fits a random forest with \(m = p/3\), reads the OOB score, computes permutation importance, and then asks the forest about a basin wetter than any it was trained on. ISLP's own lab (§8.3) does the same steps on standard datasets, and the R lab is shown in the video below.

import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.inspection import permutation_importance

rng = np.random.default_rng(42)
n = 500
P   = rng.uniform(400, 2000, n)    # annual precipitation (mm)
PET = rng.uniform(400, 1200, n)    # potential ET (mm)
F   = rng.uniform(0, 1, n)         # forest fraction (-)
Z   = rng.uniform(200, 3000, n)    # mean elevation (m): a decoy
w = 2.6                            # Fu parameter
ET = P * (1 + PET/P - (1 + (PET/P)**w)**(1/w))
Q  = P - ET - 50*F + rng.normal(0, 20, n)
X  = np.column_stack([P, PET, F, Z])
names = ["P", "PET", "forest", "elevation"]

rf = RandomForestRegressor(n_estimators=500, max_features=1/3,
                           min_samples_leaf=5, oob_score=True, random_state=0)
rf.fit(X, Q)
print(f"OOB R^2: {rf.oob_score_:.3f}")

imp = permutation_importance(rf, X, Q, n_repeats=10, random_state=0)
for name, m in sorted(zip(names, imp.importances_mean), key=lambda t: -t[1]):
    print(f"{name:>10s}: {m:.3f}")

# Extrapolation: a basin wetter than any in the training set
x_new = np.array([[3000, 800, 0.5, 1000]])
et_new = 3000 * (1 + 800/3000 - (1 + (800/3000)**w)**(1/w))
print(f"True Q: {3000 - et_new - 25:.0f} mm, RF Q: {rf.predict(x_new)[0]:.0f} mm")

The output (scikit-learn, random_state fixed):

OOB R^2: 0.896
         P: 1.219
       PET: 0.185
    forest: 0.020
 elevation: 0.020
True Q: 2212 mm, RF Q: 918 mm   (max Q in training: 1559 mm)

Three lessons in a few lines. The OOB \(R^2\) of 0.90 is a decent fit obtained with no tuning. The importance ranking recovers precipitation and PET, but the real forest effect is indistinguishable from the decoy: small physical effects drown in noise. And for \(P = 3000\) mm the forest predicts 918 mm against a true 2212 mm, below even the largest runoff in its training set. That is §6's extrapolation warning, in numbers.

The same experiment in Java, with Smile

For those of us whose models live on the JVM, as GEOframe does, Smile (Statistical Machine Intelligence and Learning Engine, by Haifeng Li) offers a random forest that can sit inside the same code base. The defaults match the book: \(m = p/3\) for regression, minimum leaf size 5, OOB metrics computed during training. Add the dependency (Smile 4 and later need Java 21 or newer):

<dependency>
  <groupId>com.github.haifengl</groupId>
  <artifactId>smile-core</artifactId>
  <version>6.3.0</version>
</dependency>
import java.util.Random;
import smile.data.DataFrame;
import smile.data.formula.Formula;
import smile.regression.RandomForest;

public class BudykoForest {

    // Fu's form of the Budyko curve: actual ET from P and PET
    static double fuET(double P, double PET, double w) {
        return P * (1 + PET / P - Math.pow(1 + Math.pow(PET / P, w), 1 / w));
    }

    public static void main(String[] args) {
        int n = 500;
        double w = 2.6;
        Random rng = new Random(42);

        double[][] data = new double[n][];
        double qMax = Double.NEGATIVE_INFINITY;
        for (int i = 0; i < n; i++) {
            double P   = 400 + 1600 * rng.nextDouble();   // annual precipitation (mm)
            double PET = 400 +  800 * rng.nextDouble();   // potential ET (mm)
            double F   = rng.nextDouble();                // forest fraction (-)
            double Z   = 200 + 2800 * rng.nextDouble();   // mean elevation (m): a decoy
            double Q   = P - fuET(P, PET, w) - 50 * F + 20 * rng.nextGaussian();
            data[i] = new double[] {P, PET, F, Z, Q};
            qMax = Math.max(qMax, Q);
        }
        String[] names = {"P", "PET", "forest", "elevation", "Q"};
        DataFrame df = DataFrame.of(data, names);

        // 500 trees; mtry defaults to p/3, nodeSize to 5
        RandomForest rf = RandomForest.fit(Formula.lhs("Q"), df,
                                           new RandomForest.Options(500));

        System.out.printf("OOB R^2: %.3f%n", rf.metrics().r2());

        // Impurity-based importance: total RSS decrease, summed over trees
        double[] imp = rf.importance();
        for (int j = 0; j < imp.length; j++) {
            System.out.printf("%10s: %.3e%n", names[j], imp[j]);
        }

        // Extrapolation: a basin wetter than any in the training set
        DataFrame wet = DataFrame.of(new double[][] {{3000, 800, 0.5, 1000, 0}}, names);
        double qTrue = 3000 - fuET(3000, 800, w) - 25;
        System.out.printf("True Q: %.0f mm, RF Q: %.0f mm (max Q in training: %.0f mm)%n",
                          qTrue, rf.predict(wet.get(0)), qMax);
    }
}

Java's random generator is not NumPy's, so the numbers will differ slightly from the Python run, but the story does not: a high OOB \(R^2\), precipitation dominating, and a wet-basin prediction capped near the training maximum. One difference is instructive. Smile's importance() is the impurity-based measure of §4, not permutation importance, so expect the decoy elevation to receive a visibly non-zero score. That is warning 3 of §4 at work.

Lab: Random Forests and Boosting (Ch. 8, 15:35), in R. A companion Lab: Decision Trees precedes it.

References and lectures

The basis of this post

For going further

  • Hastie, T., Tibshirani, R., Friedman, J. (2009). The Elements of Statistical Learning, 2nd ed. Springer. Chapter 15, Random Forests.
  • Breiman, L. (2001). Random forests. Machine Learning, 45, 5–32.
  • Strobl, C., Boulesteix, A.-L., Zeileis, A., Hothorn, T. (2007). Bias in random forest variable importance measures. BMC Bioinformatics, 8, 25.
  • Meinshausen, N. (2006). Quantile regression forests. Journal of Machine Learning Research, 7, 983–999.

Software

  • Li, H. Smile: Statistical Machine Intelligence and Learning Engine. haifengl.github.io; artifact smile-core on Maven Central.
  • Pedregosa, F. et al. (2011). Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12, 2825–2830.
Posted by at
Labels: , ,

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)