Random Forests, and Why Averaging Works

Context

A random forest (an average of many fully grown decision trees, each grown on its own random draw of the training rows and able to choose each split from a random subset of the columns) prices a listing at the mean of the prices its trees give it. One fully grown tree changes a great deal when its rows change, and the mean of many trees that change in different directions changes far less.

It is the model a pricing team reaches for when a table mixes numbers, categories and gaps and the price does not follow straight lines: a second opinion beside the linear pipeline, a read of which columns a fitted model leaned on, and, grown as a classifier, the chance that a listing sells within 30 days. Lesson 3.6.2 cross-validates a forest on the same folds as the linear pipeline, and Lesson 3.6.3 searches a tree's and a boosted model's settings on such folds.

The failures look like modeling work. Two forests built with the same seed average to exactly one of them, a hundred trees grown on the same rows gain little over one, and max_features=1 offers each split a single column where 1.0 offers all of them. A column of random numbers outranks a listing's condition in a forest's importances, and a forest that holds its own on price loses on whether a listing sells. Lesson 3.5.1 grew a tree that learned its rows by heart and changed whenever they did, and this lesson averages many such trees.

Goal

• Average two fully grown trees and account for the gain exactly, as the trees' average squared miss less their spread.
• Build a RandomForestRegressor in a pipeline, read its trees from estimators_ and their bootstrap draws from estimators_samples_, and show its prediction is the mean of its trees'.
• Choose max_features on validation rows, distinguish an int from a float and "sqrt", and explain from the trees' spread and strength when fewer columns per split help.
• Measure how the validation score settles as trees are added, fix the forest's settings, and score the test rows once.
• Reproduce one out-of-bag prediction by hand, and identify what feature_importances_ measures and what it never says.
• Compare a RandomForestClassifier with the logistic pipeline by ROC AUC and average precision, and explain its probabilities as votes.

Prerequisites

None of this is required, and the lesson stands on its own. Both are here if you want a forest built step by step, or scikit-learn's own account of one, before you start.

• YouTube - "StatQuest: Random Forests Part 1 - Building, Using and Evaluating" by StatQuest with Josh Starmer. A forest built in three steps, a bootstrapped dataset, a tree that picks each split from a random subset of the variables, and both repeated many times, then evaluated on the samples each tree left out.
• Official docs - scikit-learn 1.9 User Guide, "Ensembles", the section "Random Forests". Each tree built from a bootstrap sample with a random subset of the features considered at each split, and a classifier that averages its trees' probabilities rather than counting their votes.

Table of Contents

• Averaging Two Trees
• Bootstrap Samples
• A Random Subset of Features at Each Split
• How Many Trees
• Out-of-Bag Rows and Feature Importances
• A Forest That Classifies

Theory

Nothing to install: everything this lesson imports from scikit-learn arrived with Lesson 3.1.1's install. Every Theory block goes into one script, forests.py, saved in ~/dev/bootcamp beside listings.csv and run with python3 forests.py, each block under the one before, until the last subtitle starts a script of its own. The rows are Lesson 3.3.1's: 600 test listings held out of all 3,000, then 1,800 fit rows and 600 validation rows carved from the rest. Every setting is compared on the validation rows, and the test rows score the forest once, after its settings are fixed.

Averaging Two Trees

What. Averaging two models prices each listing at the mean of the two prices they give it. Where both miss a listing on the same side, the average misses by the mean of their two misses, and where one prices too high and the other too low, part of each miss cancels. Why. A fully grown tree learns its training rows by heart, so two trees grown on different rows disagree on many new listings, and that disagreement is what averaging turns into accuracy. Its size decides how much. When. Whenever several models of one target miss in different directions on the same listings, and above all with fully grown trees, whose misses swing with the rows they were grown on.

The two trees are Lesson 3.5.1's make_tree() grown on the first 900 fit rows and on the last 900:

import numpy as np
import pandas as pd
from sklearn.compose import ColumnTransformer
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import OneHotEncoder, OrdinalEncoder
from sklearn.tree import DecisionTreeRegressor

listings = pd.read_csv("listings.csv")
numeric = ["area_sqm", "year_built", "distance_km", "rooms", "bathrooms"]
nominal = ["neighborhood", "property_type", "energy_rating"]
levels = ["poor", "fair", "good", "excellent"]
features = numeric + nominal + ["condition"]
X_train, X_test, y_train, y_test = train_test_split(
    listings[features], listings["price"], test_size=0.2, random_state=42
)
X_fit, X_val, y_fit, y_val = train_test_split(
    X_train, y_train, test_size=0.25, random_state=42
)


def tree_prep(num_cols):
    # Lesson 3.5.1's make_prep("passthrough"), columns as a parameter
    return ColumnTransformer([
        ("num", "passthrough", num_cols),
        ("cat", OneHotEncoder(sparse_output=False), nominal),
        ("cond", OrdinalEncoder(categories=[levels]), ["condition"]),
    ])


def make_tree():
    tree = DecisionTreeRegressor(random_state=42)
    return make_pipeline(tree_prep(numeric), tree)


tree_a = make_tree().fit(X_fit.iloc[:900], y_fit.iloc[:900])
tree_b = make_tree().fit(X_fit.iloc[900:], y_fit.iloc[900:])
pred_a = tree_a.predict(X_val)
pred_b = tree_b.predict(X_val)
averaged = (pred_a + pred_b) / 2
print(round(r2_score(y_val, pred_a), 3))    # 0.648
print(round(r2_score(y_val, pred_b), 3))    # 0.607
print(round(r2_score(y_val, averaged), 3))  # 0.776

The averaged prediction scores 0.776, above both trees. Where the gain comes from, listing by listing:

prices = y_val.to_numpy()
miss_a = prices - pred_a
miss_b = prices - pred_b
apart = miss_a * miss_b < 0  # one too high, one too low
two_misses = (np.abs(miss_a) + np.abs(miss_b)) / 2
avg_miss = np.abs(prices - averaged)
print(apart.sum())           # 227
print(round(avg_miss[~apart].mean()), round(two_misses[~apart].mean()))
print(round(avg_miss[apart].mean()), round(two_misses[apart].mean()))

83563 83563
32427 75198

On the 373 listings where no tree prices too high while the other prices too low, the average misses by exactly the mean of the two misses, 83,563, and gains nothing. On the 227 missed on opposite sides, part of each miss cancels, 32,427 against 75,198. Lesson 3.1.3's split of a squared loss accounts for the whole gain at once, with a listing's two predictions in place of the prices and its own price in place of the candidate:

both = np.array([pred_a, pred_b])            # one row per tree
tree_sq = ((prices - both) ** 2).mean()
spread = np.var(both, axis=0).mean()
avg_sq = ((prices - averaged) ** 2).mean()
print(np.isclose(avg_sq, tree_sq - spread))  # True
var_val = np.var(y_val)
print(round(1 - tree_sq / var_val, 3))       # 0.628
print(round(spread / var_val, 3))            # 0.148

tree_sq is the trees' squared miss averaged over both trees and all 600 listings, and spread is how far a listing's two predictions sit from their own mean, Lesson 2.2.5's mean squared distance from the mean, averaged over the listings. The averaged prediction's squared miss is the first less the second, exactly. Divided by np.var(y_val), the validation prices' own spread, which Lesson 3.1.4's R² measures against, each term lands on R²'s scale: the trees average 0.628, their spread is worth 0.148, and the average scores the sum, 0.776. Averaging gains the spread and nothing else, so the rest of this lesson grows trees that differ.

DON'T:

• DON'T average copies of one model. Three copies of the first tree average to its own 0.648, because copies have no spread.
• DON'T expect the averaged prediction to beat the better tree on every listing. It beat both trees on 134 of the 600 listings and missed by more than the better one on 462, and still scored 0.776 against 0.648, because the better tree changes from one listing to the next.

Bootstrap Samples

What. A bootstrap sample (as many rows as the training rows hold, drawn at random with replacement, so some rows come twice or more and some never) gives each tree its own version of the rows. RandomForestRegressor, from sklearn.ensemble, grows each of its trees, 100 unless n_estimators says otherwise, fully on its own bootstrap sample, and predicts their mean. Why. Averaging gains the trees' spread, and the two halves bought spread by giving each tree half the listings. A bootstrap sample changes which listings a tree sees, and how often, while keeping as many draws as there are rows. When. Whenever a model's predictions swing with its training rows, a fully grown tree's above all.

from sklearn.ensemble import RandomForestRegressor


def make_forest(max_features=1.0, n_estimators=100, bootstrap=True):
    forest = RandomForestRegressor(
        n_estimators=n_estimators,
        max_features=max_features,
        bootstrap=bootstrap,
        random_state=42,
    )
    return make_pipeline(tree_prep(numeric), forest)


forest = make_forest().fit(X_fit, y_fit)
single = make_tree().fit(X_fit, y_fit)
print(round(forest.score(X_val, y_val), 3))  # 0.849
print(round(single.score(X_val, y_val), 3))  # 0.614
print(len(forest[-1].estimators_))           # 100
print(forest[-1].estimators_[0])

DecisionTreeRegressor(max_features=1.0, random_state=1608637542)

make_forest() with no arguments is RandomForestRegressor at its defaults, seeded with 42. One tree grown on all 1,800 fit rows scores 0.614, the forest 0.849, beside Lesson 3.3.1's linear baseline, 0.84 on the same rows. estimators_ holds the 100 fitted trees, each a DecisionTreeRegressor with its own random_state, drawn from the forest's 42, so each draws its own sample, and estimators_samples_ gives the draws back, rebuilt from each tree's random_state at every call:

drawn = forest[-1].estimators_samples_[0]
times = np.bincount(drawn, minlength=len(X_fit))
print(len(drawn), (times > 0).sum())  # 1800 1139
print(pd.Series(times).value_counts().sort_index().to_dict())

{0: 661, 1: 671, 2: 323, 3: 108, 4: 26, 5: 11}

estimators_samples_ returns each tree's draw as positions in the fit rows, tree 0's first, and np.bincount (NumPy's count of how often each whole number appears in an array) counts each position's draws, minlength keeping a count for positions never drawn. Tree 0 hit 1,139 different listings in its 1,800 draws, 671 once, 323 twice and 11 five times, and never drew the other 661. It was fitted on that draw, each listing counting as many times as it was drawn, so those 661 are listings it never saw. Each tree's predictions:

seen = forest[:-1].transform(X_val)
trees = forest[-1].estimators_
per_tree = np.array([tree.predict(seen) for tree in trees])
print(per_tree.shape)  # (100, 600)
# True: the forest's price is its trees' mean
print(np.allclose(per_tree.mean(axis=0), forest.predict(X_val)))
tree_r2 = [r2_score(y_val, row) for row in per_tree]
# 0.327 0.733
print(round(min(tree_r2), 3), round(max(tree_r2), 3))

forest[:-1], Lesson 3.2.3's slice of every step but the last, puts the validation rows in the form the trees were fitted on, and each row of per_tree is one tree's 600 predictions. Their mean is the forest's prediction, to within rounding. Alone, the trees score from 0.327 to 0.733, every one below the forest's 0.849. The first subtitle's accounting, for this forest and for one whose trees all see the same rows:

def forest_report(fitted):
    # the forest's R², its trees' average R² and their spread
    seen = fitted[:-1].transform(X_val)
    trees = fitted[-1].estimators_
    per_tree = np.array([tree.predict(seen) for tree in trees])
    tree_sq = ((prices - per_tree) ** 2).mean()
    spread = np.var(per_tree, axis=0).mean()
    three = [fitted.score(X_val, y_val), 1 - tree_sq / var_val, spread / var_val]
    return np.round(three, 3)


same_rows = make_forest(bootstrap=False).fit(X_fit, y_fit)
print(forest_report(forest))     # [0.849 0.597 0.253]
print(forest_report(same_rows))  # [0.663 0.581 0.082]
seen_fit = same_rows[:-1].transform(X_fit)
fit_preds = [tree.predict(seen_fit) for tree in same_rows[-1].estimators_]
# True: every tree returns every fit row's own price
print(np.all(np.array(fit_preds) == y_fit.to_numpy()))

The forest's R² is its trees' average R² plus their spread, up to rounding in the third decimal. With bootstrap=False every tree is fitted on all 1,800 rows with every column and returns every fit row's own price, so the trees agree on every row they were fitted on. They still differ on new listings, because each tries the columns in its own random order, and scikit-learn settles two splits that lower the squared error equally by that order. Their spread is worth 0.082, a third of the bootstrapped trees', each tree is about as good, 0.581 against 0.597, and the forest scores 0.663.

DON'T:

• DON'T switch bootstrap off to give every tree every row. Each tree barely changed, 0.581 against 0.597 on average, and the forest fell from 0.849 to 0.663, because trees fitted on the same rows spread a third as far.
• DON'T judge a forest by one of its trees. The 100 trees scored from 0.327 to 0.733 on the validation rows, every one below the forest's 0.849.

A Random Subset of Features at Each Split

What. max_features sets how many of the 24 columns the preparation hands the forest each split may choose among, drawn afresh at every split. A float is a share of the columns, an int a count of them, and "sqrt" the square root of their number, rounded down. RandomForestRegressor defaults to 1.0, every column at every split, and RandomForestClassifier to "sqrt". Why. With every column on offer, a tree's first split is the one that lowers the squared error most on its draw, and on these listings that is nearly always a question about area, so the trees open alike. Offered fewer columns, a split must sometimes do without area, and the trees differ more, at a cost to each tree. When. For every forest, as a setting chosen on validation rows. Lesson 3.6.3 searches settings like it together, on folds.

Four settings, from every column to a single one, are scored on the validation rows by a rule written down before the loop: the highest validation R² to three decimals.

names = forest[:-1].get_feature_names_out()
val_r2 = {}
for setting in [1.0, 0.5, "sqrt", 1]:
    fitted = make_forest(max_features=setting).fit(X_fit, y_fit)
    trees = fitted[-1].estimators_
    # tree_ is a tree's fitted structure, and tree_.feature holds the
    # position of the column each node's split asks about, the root first
    roots = [names[tree.tree_.feature[0]] for tree in trees]
    report = forest_report(fitted)
    val_r2[setting] = report[0]
    print(setting, trees[0].max_features_, roots.count("num__area_sqm"), report)
print(max(val_r2, key=val_r2.get))  # 0.5

1.0 24 99 [0.849 0.597 0.253]
0.5 12 48 [0.859 0.593 0.267]
sqrt 4 17 [0.85  0.526 0.324]
1 1 5 [0.801 0.378 0.424]

max_features_, set by each tree's fit, is how many columns each of its splits chose among, and each line holds the setting, that number, how many of the 100 trees open on area, and forest_report's three numbers.

Area is among 12 columns drawn from 24 half the time, among 4 a sixth of the time and alone a 24th of the time, and the trees open on it about that often, 48, 17 and 5 times in 100, against 99 with every column. As fewer columns are offered, the trees spread wider, from 0.253 to 0.424, and each gets worse alone, from 0.597 to 0.378. The forest gains where the added spread outweighs what each tree loses: at 0.5 the spread grew by 0.014 while each tree lost 0.004, and the forest rose to 0.859. At "sqrt" both moved by 0.071 and the forest stayed at 0.85, and at a single column each tree lost 0.219 for 0.171 more spread, and it fell to 0.801. The runs show that trade, not which listings the trees disagree on or why. The rule picks 0.5.

DON'T:

• DON'T write max_features=1 meaning every column. An int is a count, so 1 offered each split a single column and scored 0.801, where the float 1.0 offered all 24 and scored 0.849.
• DON'T carry "sqrt" over from the classifier as the regressor's default. RandomForestRegressor defaults to 1.0, and on these rows "sqrt", 4 columns per split, scored 0.85, below 0.5's 0.859.

How Many Trees

What. n_estimators sets how many trees the forest grows, 100 unless set. Each tree added moves the forest's mean prediction by the tree's own distance from the old mean, divided by the new count of trees, so the more trees a forest holds, the less a tree at a given distance from the mean moves it. Why. Every tree is fitted and every tree predicts every listing, so the work grows with the count while the score stops rising once the mean has settled. The count is set where more trees stop moving the validation score. When. After the other settings are chosen, and before the final fit.

One forest of 300 trees at the chosen max_features, beside the 100-tree forest:

many = make_forest(max_features=0.5, n_estimators=300).fit(X_fit, y_fit)
hundred = make_forest(max_features=0.5).fit(X_fit, y_fit)
seen = many[:-1].transform(X_val)
per_tree = np.array([tree.predict(seen) for tree in many[-1].estimators_])
# True: the first 100 of the 300 trees are the 100-tree forest
print(np.allclose(per_tree[:100].mean(axis=0), hundred.predict(X_val)))
for n in [1, 10, 50, 100, 150, 200, 300]:
    print(n, round(r2_score(y_val, per_tree[:n].mean(axis=0)), 3))

1 0.447
10 0.82
50 0.857
100 0.859
150 0.862
200 0.862
300 0.861

The first 100 of the 300 trees are the 100-tree forest's, so the mean of the first n is the forest of n trees. From 50 trees on the score stays between 0.857 and 0.862. How far each stretch of trees moves an average listing's prediction, and the score on the fit rows:

for start, end in [(10, 30), (30, 100), (100, 300)]:
    early = per_tree[:start].mean(axis=0)
    late = per_tree[:end].mean(axis=0)
    print(start, end, round(np.abs(late - early).mean()))
print(round(hundred.score(X_fit, y_fit), 3))
print(round(many.score(X_fit, y_fit), 3))

10 30 17716
30 100 10371
100 300 5812
0.969
0.97

Each stretch moves the mean less than the one before, though it holds more trees, and on the fit rows the score moves only from 0.969 to 0.97: the extra trees settle the mean and barely change how closely the forest follows its rows, where Lesson 3.5.3's boosting will differ. Past 100 trees the validation score moves by at most 0.003 for three times the trees, so the forest keeps 100. Its settings are fixed, the validation rows rejoin the training rows, and the test rows score it once:

final = make_forest(max_features=0.5).fit(X_train, y_train)
test_pred = final.predict(X_test)
test_r2 = final.score(X_test, y_test)
print(np.round(test_pred[:3]))  # [196600. 457350. 551150.]
print(round(test_r2, 3))        # 0.775

The forest scores 0.775 on the 600 test listings, where Lesson 3.2.3's linear pipeline scored 0.766 and Lesson 3.5.1's fully grown tree 0.587, its chosen tree 0.673. Its 0.859 on the validation rows is not set beside this 0.775, Lesson 3.3.1's rule, and its 0.009 over the linear pipeline is one draw of 600 listings, the kind of gap Lesson 3.7.1 tests.

DON'T:

• DON'T expect more trees to keep raising the score. From 150 trees to 300 the validation R² went 0.862, 0.862 and 0.861.
• DON'T cut a forest to a handful of trees unmeasured. Ten trees scored 0.82 on the validation rows, against 0.857 for fifty.

Out-of-Bag Rows and Feature Importances

What. A tree's out-of-bag rows (the training rows its bootstrap sample never drew) are rows it was never fitted on. With oob_score=True the forest predicts each training row with only the trees that never drew it, and keeps those predictions in oob_prediction_ and their R² in oob_score_. A feature importance (each tree's reduction of the squared error on the rows it drew, credited to the column each split asked about, turned into shares and averaged over the trees) is kept in feature_importances_, the shares summing to 1. Why. Every training row is new to about a third of the trees, so the out-of-bag score is a held-out score that sets no rows aside. The importances say which columns the splits used on the rows the trees drew, which is not what moves a price. When. The out-of-bag score when no rows can be spared for validation, set beside nothing but another out-of-bag score. The importances when reading what a fitted forest leaned on, never as an account of a price.

The chosen forest grown on the fit rows again, its out-of-bag predictions kept:

bagged = RandomForestRegressor(
    max_features=0.5, oob_score=True, random_state=42
)
oob_forest = make_pipeline(tree_prep(numeric), bagged).fit(X_fit, y_fit)
print(round(bagged.oob_score_, 3))  # 0.771
# 0.859
print(round(oob_forest.score(X_val, y_val), 3))
out_of_bag = np.zeros(len(X_fit), dtype=int)
for drawn in bagged.estimators_samples_:
    out_of_bag += np.bincount(drawn, minlength=len(X_fit)) == 0
print(out_of_bag.min(), out_of_bag.max(), round(out_of_bag.mean(), 1))

21 53 36.8

These are the 100 trees hundred holds, so the validation score is its 0.859: keeping the out-of-bag predictions changes nothing in the fit. Each fit row is out of bag for 21 to 53 of the trees, 36.8 on average. The first fit row, by hand:

first_row = oob_forest[:-1].transform(X_fit.iloc[[0]])
never_drew = [
    tree.predict(first_row)[0]
    for tree, drawn in zip(bagged.estimators_, bagged.estimators_samples_)
    if 0 not in drawn
]
# 32 371656
print(len(never_drew), round(np.mean(never_drew)))
print(round(bagged.oob_prediction_[0]))  # 371656

The 32 trees that never drew the first fit row average 371,656, its out-of-bag prediction. oob_score_ is the R² of all 1,800 such predictions, 0.771, where the same forest scores 0.859 on the validation rows: two honest held-out scores on two sets of listings, which differ as two draws do, and by a little more, because each out-of-bag prediction averages only the trees that never drew its row, about a third of them. The final forest's importances:

names = final[:-1].get_feature_names_out()
shares = pd.Series(final[-1].feature_importances_, index=names)
print(shares.sum())  # 1.0
print(shares.sort_values(ascending=False).head(5).round(3))

num__area_sqm                 0.417
num__rooms                    0.162
num__distance_km              0.102
cat__neighborhood_Eastgate    0.074
num__year_built               0.048
dtype: float64

On average over the trees, splits on area removed 0.417 of the squared error a tree's splits removed, splits on rooms 0.162 and on distance 0.102. Two runs show what such a share does not say: a column of random numbers, drawn with Lesson 2.2.1's default_rng(42), joins the training rows, and the forest that offers every column at every split is grown beside the final one:

rng = np.random.default_rng(42)
noisy = X_train.assign(noise=rng.normal(0, 1, len(X_train)))
with_noise = make_pipeline(
    tree_prep(numeric + ["noise"]),
    RandomForestRegressor(max_features=0.5, random_state=42),
).fit(noisy, y_train)
noise_names = with_noise[:-1].get_feature_names_out()
added = pd.Series(with_noise[-1].feature_importances_, index=noise_names)
place = list(added.sort_values(ascending=False).index).index("num__noise")
print(place + 1, len(added))               # 9 25
print(round(added["num__noise"], 3))       # 0.027
print(round(added["cond__condition"], 3))  # 0.022
every_column = make_forest().fit(X_train, y_train)
all_shares = pd.Series(every_column[-1].feature_importances_, index=names)
print(round(all_shares["num__rooms"], 3))  # 0.029

The random numbers cannot move a price, and they took 0.027, ninth of 25 columns, ahead of condition's 0.022. Rooms took 0.162 in the final forest and 0.029 in the forest that offers every column, on the same rows and the same prices. An importance is a fact about the splits on the rows the trees drew: random numbers earned a share by being split on, and rooms' share moved more than fivefold with a setting that only changes which columns a split may choose among. Neither share says what one more room would do to a listing's price.

DON'T:

• DON'T read an importance as what a column does to a price. A column of random numbers took 0.027, ahead of condition's 0.022, and rooms took 0.162 with half the columns per split and 0.029 with all of them.
• DON'T set an out-of-bag score beside a validation score. The same forest scored 0.771 out of bag and 0.859 on the validation rows, two different sets of listings, each out-of-bag prediction made by about a third of the trees.
• DON'T read oob_score_ from a handful of trees. With 5 trees, 188 of the 1,800 fit rows were drawn by all five, scikit-learn warned Some inputs do not have OOB scores. This probably means too few trees were used to compute any reliable OOB estimates., those rows got an out-of-bag prediction of 0, and the score came out 0.127.

A Forest That Classifies

What. RandomForestClassifier, from sklearn.ensemble, grows classification trees on bootstrap samples, and its predict_proba averages its trees' probabilities. On these listings each fully grown tree's probability of sold is 0 or 1, so the forest's is the share of its trees that vote sold. Its max_features defaults to "sqrt". Why. Averaging works the same way on a label, and whether the forest beats the logistic pipeline is a question for the validation rows, read with Section 3.4's tools. A forest's showing on one target says nothing about the next. When. For a category target, with the forest set against the logistic pipeline on the same rows before either is trusted.

The subtitle starts sold_forest.py, on the sold label's own split, Section 3.4's, stratified on the label, with a validation carve stratified the same way:

import numpy as np
import pandas as pd
from sklearn.compose import ColumnTransformer
from sklearn.ensemble import RandomForestClassifier
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import average_precision_score, roc_auc_score
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import OneHotEncoder, OrdinalEncoder, StandardScaler

listings = pd.read_csv("listings.csv")
numeric = ["area_sqm", "year_built", "distance_km", "rooms", "bathrooms"]
numeric = numeric + ["price"]  # known when the listing goes up
nominal = ["neighborhood", "property_type", "energy_rating"]
levels = ["poor", "fair", "good", "excellent"]
X = listings[numeric + nominal + ["condition"]]
y = listings["sold_in_30_days"]
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y
)
X_fit, X_val, y_fit, y_val = train_test_split(
    X_train, y_train, test_size=0.25, random_state=42, stratify=y_train
)


def prep(num_step):
    return ColumnTransformer([
        ("num", num_step, numeric),
        ("cat", OneHotEncoder(sparse_output=False), nominal),
        ("cond", OrdinalEncoder(categories=[levels]), ["condition"]),
    ])


def logistic():
    fill = SimpleImputer(strategy="median", add_indicator=True)
    scaled = make_pipeline(fill, StandardScaler())
    return make_pipeline(prep(scaled), LogisticRegression())


def voting_forest():
    forest = RandomForestClassifier(random_state=42)
    return make_pipeline(prep("passthrough"), forest)


def two_scores(make, X_a, y_a, X_b, y_b):
    sold = make().fit(X_a, y_a).predict_proba(X_b)[:, 1]
    return [roc_auc_score(y_b, sold), average_precision_score(y_b, sold)]


for name, make in [("logistic", logistic), ("forest", voting_forest)]:
    scores = two_scores(make, X_fit, y_fit, X_val, y_val)
    print(name, np.round(scores, 3))

logistic [0.647 0.289]
forest [0.602 0.251]

two_scores fits a pipeline on one set of rows and scores its probability of sold, column 1 of Lesson 3.4.1's predict_proba, on another, by Lesson 3.4.3's ROC AUC and average precision. On the validation rows the logistic pipeline leads on both. Whether the listings a split held out decided the lead, Lesson 3.3.1 counted across 20 splits of the training rows, each stratified here on the label:

logistic_ahead = np.zeros(2, dtype=int)
for seed in range(20):
    X_a, X_b, y_a, y_b = train_test_split(
        X_train, y_train, test_size=0.25, random_state=seed, stratify=y_train
    )
    lin = two_scores(logistic, X_a, y_a, X_b, y_b)
    rf = two_scores(voting_forest, X_a, y_a, X_b, y_b)
    logistic_ahead += np.array(lin) > np.array(rf)
print(logistic_ahead)  # [20 20]

The logistic pipeline led on both measures on all 20 splits, so which listings a split held out did not decide the lead. On this label the forest at its defaults is the weaker model. Each pipeline is refitted on all 2,400 training rows and scores the test rows once:

for name, make in [("logistic", logistic), ("forest", voting_forest)]:
    scores = two_scores(make, X_train, y_train, X_test, y_test)
    print(name, np.round(scores, 3))

logistic [0.642 0.306]
forest [0.604 0.278]

The logistic pipeline's 0.642 and 0.306 are the scores Section 3.4 gave it on these 600 test listings, and the forest stays behind on both. Where the forest's probabilities come from:

forest = voting_forest().fit(X_train, y_train)
seen = forest[:-1].transform(X_test)
votes = np.array([tree.predict(seen) for tree in forest[-1].estimators_])
print(np.unique(votes))          # [0. 1.]
share = votes.mean(axis=0)       # the share of trees voting sold
sold = forest.predict_proba(X_test)[:, 1]
print(np.allclose(share, sold))  # True
print(np.round(sold[:3], 2))     # [0.07 0.33 0.17]
train_auc = roc_auc_score(y_train, forest.predict_proba(X_train)[:, 1])
print(round(train_auc, 3))       # 1.0

Inside the forest each tree predicts a class by its position in classes_, [False True], so every vote is 0.0 or 1.0, and the mean of the 100 votes is predict_proba's column for True: 7 of the 100 trees vote sold for the first test listing. On the training rows the probabilities rank every sold listing above every unsold one, an AUC of 1.0, Lesson 3.5.1's fully grown trees learning their rows by heart, here on a label.

DON'T:

• DON'T carry a forest's showing on one target over to another. On sold_in_30_days the logistic pipeline led on ROC AUC and average precision on all 20 splits of the training rows and on the test rows, 0.642 and 0.306 against 0.604 and 0.278.
• DON'T read the forest's training AUC as skill. It is 1.0 on the rows its trees were grown on and 0.604 on the test rows.

Lab Exercises

Lab Exercise 1: Odd and Even Rows, and Two Ways for Trees to Differ

In differ.py in ~/dev/bootcamp, build the Theory's split, carve and tree preparation. Grow one fully grown tree on the fit rows at even positions, X_fit.iloc[::2], and one on those at odd positions, X_fit.iloc[1::2]. Print each tree's validation R² and the averaged prediction's, the number of validation listings the two miss on opposite sides, the check that the averaged prediction's squared miss is the trees' average squared miss less their spread, and the spread in R² units.

Then grow four forests of 100 trees on the fit rows, bootstrap at True and at False, each with max_features at 1.0 and at 0.5, and print for each the two settings and the Theory's three numbers: the forest's validation R², its trees' average R² and their spread. End with a comment saying which of the two ways to make trees differ did more on these rows.

Lab Exercise 2: Out of Bag at Every Size, and Importances Without Rooms

In oob_sizes.py, on the Theory's split and carve, grow forests of 25, 50, 100, 200 and 400 trees on the fit rows with max_features=0.5 and oob_score=True, and print each count beside its out-of-bag R² and its validation R². For the 400-tree forest, print the fewest, the average and the most trees any fit row is out of bag for.

Then grow the Theory's final forest on all 2,400 training rows twice, once on the Theory's columns and once with rooms left out of the numeric ones, and print each forest's five largest importances. Say in a comment where rooms' share went, and what that says about reading an importance.

import numpy as np
import pandas as pd
from sklearn.compose import ColumnTransformer
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import OneHotEncoder, OrdinalEncoder
from sklearn.tree import DecisionTreeRegressor

listings = pd.read_csv("listings.csv")
numeric = ["area_sqm", "year_built", "distance_km", "rooms", "bathrooms"]
nominal = ["neighborhood", "property_type", "energy_rating"]
levels = ["poor", "fair", "good", "excellent"]
features = numeric + nominal + ["condition"]
X_train, X_test, y_train, y_test = train_test_split(
    listings[features], listings["price"], test_size=0.2, random_state=42
)
X_fit, X_val, y_fit, y_val = train_test_split(
    X_train, y_train, test_size=0.25, random_state=42
)
prices = y_val.to_numpy()
var_val = np.var(y_val)


def tree_prep(num_cols):
    return ColumnTransformer([
        ("num", "passthrough", num_cols),
        ("cat", OneHotEncoder(sparse_output=False), nominal),
        ("cond", OrdinalEncoder(categories=[levels]), ["condition"]),
    ])


def make_tree():
    return make_pipeline(tree_prep(numeric), DecisionTreeRegressor(random_state=42))


even = make_tree().fit(X_fit.iloc[::2], y_fit.iloc[::2])
odd = make_tree().fit(X_fit.iloc[1::2], y_fit.iloc[1::2])
both = np.array([even.predict(X_val), odd.predict(X_val)])
averaged = both.mean(axis=0)
for pred in [both[0], both[1], averaged]:
    print(round(r2_score(y_val, pred), 3))
misses = prices - both
tree_sq = (misses ** 2).mean()
spread = np.var(both, axis=0).mean()
avg_sq = ((prices - averaged) ** 2).mean()
print((misses[0] * misses[1] < 0).sum())     # 196
print(np.isclose(avg_sq, tree_sq - spread))  # True
print(round(spread / var_val, 3))            # 0.107

0.657
0.644
0.758

def report(fitted):
    seen = fitted[:-1].transform(X_val)
    trees = fitted[-1].estimators_
    per_tree = np.array([tree.predict(seen) for tree in trees])
    tree_sq = ((prices - per_tree) ** 2).mean()
    spread = np.var(per_tree, axis=0).mean()
    three = [fitted.score(X_val, y_val), 1 - tree_sq / var_val, spread / var_val]
    return np.round(three, 3)


for bootstrap in [True, False]:
    for max_features in [1.0, 0.5]:
        forest = RandomForestRegressor(
            max_features=max_features, bootstrap=bootstrap, random_state=42
        )
        fitted = make_pipeline(tree_prep(numeric), forest).fit(X_fit, y_fit)
        print(bootstrap, max_features, report(fitted))
# Either way alone does most of the work: half the columns per split
# with every row, 0.848, nearly matches bootstrap samples with every
# column, 0.849, and the two together score 0.859. Without either,
# the trees spread 0.082 and the forest scores 0.663.

True 1.0 [0.849 0.597 0.253]
True 0.5 [0.859 0.593 0.267]
False 1.0 [0.663 0.581 0.082]
False 0.5 [0.848 0.623 0.224]

def half_forest(num_cols, n_estimators=100, oob_score=False):
    prep = ColumnTransformer([
        ("num", "passthrough", num_cols),
        ("cat", OneHotEncoder(sparse_output=False), nominal),
        ("cond", OrdinalEncoder(categories=[levels]), ["condition"]),
    ])
    forest = RandomForestRegressor(
        n_estimators=n_estimators,
        max_features=0.5,
        oob_score=oob_score,
        random_state=42,
    )
    return make_pipeline(prep, forest)


for n in [25, 50, 100, 200, 400]:
    fitted = half_forest(numeric, n, oob_score=True).fit(X_fit, y_fit)
    oob = fitted[-1].oob_score_
    print(n, round(oob, 3), round(fitted.score(X_val, y_val), 3))

out_of_bag = np.zeros(len(X_fit), dtype=int)
# the draws of the last forest fitted, the 400-tree one
for drawn in fitted[-1].estimators_samples_:
    out_of_bag += np.bincount(drawn, minlength=len(X_fit)) == 0
print(out_of_bag.min(), round(out_of_bag.mean(), 1), out_of_bag.max())

25 0.749 0.844
50 0.761 0.857
100 0.771 0.859
200 0.777 0.862
400 0.778 0.862
112 147.2 180

no_rooms = [col for col in numeric if col != "rooms"]
for num_cols in [numeric, no_rooms]:
    fitted = half_forest(num_cols).fit(X_train, y_train)
    names = fitted[:-1].get_feature_names_out()
    shares = pd.Series(fitted[-1].feature_importances_, index=names)
    print(shares.sort_values(ascending=False).head(5).round(3))
# Without rooms, area's share rose from 0.417 to 0.473, and the
# apartment column and bathrooms joined the five largest. A share is
# what the splits used, and it moves with the columns on offer.

num__area_sqm                 0.417
num__rooms                    0.162
num__distance_km              0.102
cat__neighborhood_Eastgate    0.074
num__year_built               0.048
dtype: float64
num__area_sqm                   0.473
num__distance_km                0.099
cat__property_type_apartment    0.079
cat__neighborhood_Eastgate      0.076
num__bathrooms                  0.056
dtype: float64

Theory Assessment

Q1: Two Trees, Three Numbers

Lesson 3.5.1's two fully grown trees scored 0.648 and 0.607 on the validation rows and their average 0.776, though the average missed by more than the better tree on 462 of the 600 listings. Explain where the gain came from, using the 227 listings the two trees missed on opposite sides and the 373 others, and say what the accounting, avg_sq equal to tree_sq - spread, adds that the two groups of listings cannot. Then say what three copies of the first tree score averaged, and why.

Hint: Write a listing's averaged miss from its two misses, once with signs that agree and once with signs that differ. For the copies, ask which term of the accounting is zero.

Q2: Fewer Columns, a Better Forest

With every column on offer, 99 of the 100 trees opened on area. Explain why 48 did at max_features=0.5, 17 at "sqrt" and 5 at 1, and use the table's average tree R² and spread to explain why 0.5 scored 0.859, "sqrt" 0.85 and 1 0.801. Then say what bootstrap=False with every column did to the spread and to the forest, and why its trees differed at all.

Hint: Ask how often one column is among 12, 4 or 1 drawn from 24. For the scores, the forest's R² is its trees' average R² plus their spread.

Q3: Settled, Out of Bag, and Split On

Explain why the 200 trees added after the first 100 moved an average validation prediction by 5,812 while the 20 added after the first 10 moved it by 17,716, and why the forest keeps 100. Then explain how the first fit row's out-of-bag prediction, 371,656, was made, why the same forest scores 0.771 out of bag and 0.859 on the validation rows, and why a column of random numbers took 0.027 of the importance. Finish with what the classifier's 0.07 for the first test listing counts.

Hint: Each added tree moves the mean by its distance from the mean over the new count. For the importance, ask what the forest credits to a column, and on which rows.

Q4: Two Shortcuts to a Forest

Two teammates want a forest's gain without the work. One of them writes two_ways.py, which opens with every line of Solution 1's differ.py above def make_tree(): and goes on:

# teammate A: a hundred trees, a hundred seeds, the same fit rows
preds = []
for seed in range(100):
    tree = DecisionTreeRegressor(random_state=seed)
    fitted = make_pipeline(tree_prep(numeric), tree).fit(X_fit, y_fit)
    preds.append(fitted.predict(X_val))
print(round(r2_score(y_val, np.mean(preds, axis=0)), 3))

# teammate B: two forests, averaged
first = make_pipeline(tree_prep(numeric), RandomForestRegressor(random_state=42))
second = make_pipeline(tree_prep(numeric), RandomForestRegressor(random_state=42))
first.fit(X_fit, y_fit)
second.fit(X_fit, y_fit)
pair = (first.predict(X_val) + second.predict(X_val)) / 2
print(round(r2_score(y_val, pair), 3))

Teammate A expects a hundred seeds to give a hundred different trees and a forest's gain. Teammate B expects two forests averaged to beat one. Read the script without running it, say what each print shows, using the numbers this lesson printed, and say where each expectation goes wrong.

Hint: Ask what a seed can still decide in a tree grown on every fit row with every column, and which of the Theory's forests teammate A rebuilt. For B, ask what random_state=42 decides in a forest.

0.659
0.849

same = np.array_equal(first.predict(X_val), second.predict(X_val))
print(same)  # True

Assignment

Assignment: What the Price Leans On, for the Sellers' Page

The listings site wants to tell each seller which details the price model leans on most, and a teammate proposes the forest's feature_importances_. The reviewer wants that ranking set against one measured on listings the forest never saw. permutation_importance, from sklearn.inspection, gives one: permutation_importance(model, X, y, n_repeats=10, random_state=42) shuffles one column's values among the rows of X, scores the fitted model on the shuffled rows, does that 10 times for every column, and returns a result whose importances_mean holds, for each column of X in order, the average fall of the model's score, R² for a regressor. Handed a pipeline and its raw columns, it shuffles the raw columns, so a neighborhood is shuffled as one column, not seven.

Write lean_on.py. Before the split, add a column of random numbers, noise, to every listing with np.random.default_rng(42), as Lesson 3.3.2 drew its noise columns. Split and carve as the Theory did, grow the Theory's chosen forest on the fit rows with noise among its numeric columns, and print its validation R² beside the same forest's without noise. Then rank the ten raw columns three ways in one table: the forest's importances summed back to the raw column each encoded column came from, the permutation importance on the validation rows, and the permutation importance on the fit rows, and print where noise stands in each. End with comments saying which ranking the sellers' page should use, how the three treat noise, and what even the best of them does not say.

Hints:

• An encoded column's name is its group's name, two underscores and the rest, and the rest is the raw column's name or starts with it and an underscore, as cat__neighborhood_Eastgate starts with neighborhood_. A Series' groupby also takes a function, which it calls on every label, and groups the values by what the function returns.
• Shuffling on the fit rows measures what the forest leans on for rows its trees drew, and on the validation rows what it leans on for new ones, so compare noise across the two.

import numpy as np
import pandas as pd
from sklearn.compose import ColumnTransformer
from sklearn.ensemble import RandomForestRegressor
from sklearn.inspection import permutation_importance
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import OneHotEncoder, OrdinalEncoder

listings = pd.read_csv("listings.csv")
rng = np.random.default_rng(42)
# random numbers read nothing from the listings, so drawing them
# before the split hands the fit nothing about a held-out row
table = listings.assign(noise=rng.normal(0, 1, len(listings)))
numeric = ["area_sqm", "year_built", "distance_km", "rooms", "bathrooms"]
nominal = ["neighborhood", "property_type", "energy_rating"]
levels = ["poor", "fair", "good", "excellent"]
columns = numeric + ["noise"] + nominal + ["condition"]
train, test = train_test_split(table, test_size=0.2, random_state=42)
fit_rows, val_rows = train_test_split(train, test_size=0.25, random_state=42)
X_fit, y_fit = fit_rows[columns], fit_rows["price"]
X_val, y_val = val_rows[columns], val_rows["price"]


def chosen_forest(num_cols):
    prep = ColumnTransformer([
        ("num", "passthrough", num_cols),
        ("cat", OneHotEncoder(sparse_output=False), nominal),
        ("cond", OrdinalEncoder(categories=[levels]), ["condition"]),
    ])
    forest = RandomForestRegressor(max_features=0.5, random_state=42)
    return make_pipeline(prep, forest)


model = chosen_forest(numeric + ["noise"]).fit(X_fit, y_fit)
plain = chosen_forest(numeric).fit(X_fit, y_fit)
print(round(model.score(X_val, y_val), 3))  # 0.857
print(round(plain.score(X_val, y_val), 3))  # 0.859

def raw_column(name):
    encoded = name.split("__")[1]  # cat__neighborhood_Eastgate
    for col in columns:
        if encoded == col or encoded.startswith(col + "_"):
            return col


names = model[:-1].get_feature_names_out()
impurity = pd.Series(model[-1].feature_importances_, index=names)
on_val = permutation_importance(
    model, X_val, y_val, n_repeats=10, random_state=42
)
on_fit = permutation_importance(
    model, X_fit, y_fit, n_repeats=10, random_state=42
)
ranking = pd.DataFrame({
    "splits": impurity.groupby(raw_column).sum(),
    "shuffled_val": pd.Series(on_val.importances_mean, index=columns),
    "shuffled_fit": pd.Series(on_fit.importances_mean, index=columns),
})
print(ranking.sort_values("shuffled_val", ascending=False).round(3))
print(ranking.rank(ascending=False).loc["noise"].to_dict())
# The page uses shuffled_val, measured on 600 listings the forest never
# saw: area first, then neighborhood, distance and rooms. Shuffling
# noise there cost 0.001, last of ten, where the splits gave it 0.031,
# seventh, and shuffling it on the fit rows cost 0.027, seventh again:
# the forest used the random numbers to fit the rows its trees drew,
# and barely leans on them for new ones. Even shuffled_val says how far
# the forest's score on these listings falls when a column's values
# stop matching them, never what one more square meter or a different
# neighborhood would do to a listing's price.

               splits  shuffled_val  shuffled_fit
area_sqm        0.392         0.747         0.769
neighborhood    0.158         0.291         0.324
distance_km     0.089         0.078         0.140
rooms           0.164         0.076         0.128
condition       0.022         0.020         0.034
year_built      0.048         0.015         0.052
property_type   0.046         0.007         0.022
energy_rating   0.023         0.003         0.026
bathrooms       0.026         0.002         0.015
noise           0.031         0.001         0.027
{'splits': 7.0, 'shuffled_val': 10.0, 'shuffled_fit': 7.0}

Additional Resources

Optional, and nothing later depends on either. Worth opening for what else a forest can do, or for the two kinds of importance compared on another table.

• YouTube - "StatQuest: Random Forests Part 2: Missing data and clustering" by StatQuest with Josh Starmer. How a random forest can fill in missing values from how often two samples end in the same leaf, and how the same forest can be used to cluster samples.
• Official docs - scikit-learn 1.9 example, "Permutation Importance vs Random Forest Feature Importance (MDI)". The two importances side by side on a classifier of the Titanic's passengers, with two random columns added, and the two reasons the impurity-based one ranks a random column high: it favors columns with many distinct values, and it is computed on the training rows.

Summary

A forest's whole gain is one number, the spread of its trees, and every setting in this lesson was a way of buying spread at the least cost to each tree. Two trees grown on two halves of the rows showed the gain, bootstrap samples made it cheap, and holding columns back made trees disagree where they had opened alike. The accounting never changed, which is why copies of one tree, or a hundred trees grown on the same rows, gain little or nothing.

The numbers a forest offers about itself need the care any score needs. Its out-of-bag score is honest and costs no rows, and it is still a score of other listings than the validation rows. Its importances count what the splits used on the rows the trees drew, so random numbers earn a place among them, and a setting that only changes which columns are on offer moves a column's share more than fivefold.

The forest also lost a question, falling behind the logistic pipeline on whether a listing sells on every one of 20 splits, so no model is the answer until the rows have been asked.

Cheat Sheet

Part A - Reference Table

Part B - Compact Code Block

from sklearn.ensemble import RandomForestRegressor

forest = RandomForestRegressor(max_features=0.5, random_state=42)
pipe = make_pipeline(tree_prep(numeric), forest).fit(X_fit, y_fit)
pipe.score(X_val, y_val)           # choose settings here
seen = pipe[:-1].transform(X_val)  # the table the trees see
per_tree = np.array([t.predict(seen) for t in forest.estimators_])
per_tree.mean(axis=0)              # equals pipe.predict(X_val)
forest.estimators_samples_[0]      # tree 0's bootstrap draws
forest.feature_importances_        # split shares, summing to 1

bagged = RandomForestRegressor(oob_score=True, random_state=42)
make_pipeline(tree_prep(numeric), bagged).fit(X_fit, y_fit)
bagged.oob_score_                  # R² of out-of-bag predictions