Trees are the other family worth knowing. They find interactions without being told, need no scaling, handle mixed types, and on ordinary tabular data they are usually what wins.
One tree
A decision tree splits the data repeatedly, choosing at each step the feature and threshold that best separate the labels.
import numpy as np
import pandas as pd
from sklearn.tree import DecisionTreeClassifier, export_text
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
rng = np.random.default_rng(3)
n = 2000
tenure = rng.integers(1, 60, n)
tickets = rng.poisson(1.6, n)
spend = rng.gamma(3, 25, n).round(2)
# An INTERACTION: short tenure matters much more for low spenders.
risk = np.where(spend < 60, 0.8 - 0.015 * tenure, 0.25) + 0.08 * tickets
churned = (rng.random(n) < np.clip(risk, 0.02, 0.95)).astype(int)
X = pd.DataFrame({'tenure': tenure, 'tickets': tickets, 'spend': spend})
y = churned
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)
tree = DecisionTreeClassifier(max_depth=3, random_state=0).fit(X_tr, y_tr)
print(f'accuracy {accuracy_score(y_te, tree.predict(X_te)):.3f}')
print()
print(export_text(tree, feature_names=list(X.columns)))
Open this example in the editor and change it. Nothing you type is saved.
That printed tree is the whole model. You can read it, argue with it, and implement it as a chain of if statements. No other model family gives you that.
A tree needs no scaling and no encoding of ordinal features, because a split is a comparison rather than an arithmetic operation.
The problem is that one tree, left alone, memorises:
import numpy as np
import pandas as pd
from sklearn.tree import DecisionTreeClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
rng = np.random.default_rng(5)
n = 1500
X = rng.normal(size=(n, 6))
y = (X[:, 0] + 0.7 * X[:, 1] + rng.normal(0, 1.1, n) > 0).astype(int)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)
rows = []
for depth in [2, 3, 5, 8, 12, None]:
tree = DecisionTreeClassifier(max_depth=depth, random_state=0).fit(X_tr, y_tr)
rows.append({
'max_depth': str(depth),
'leaves': tree.get_n_leaves(),
'train': round(accuracy_score(y_tr, tree.predict(X_tr)), 3),
'test': round(accuracy_score(y_te, tree.predict(X_te)), 3),
})
table = pd.DataFrame(rows)
table['gap'] = (table['train'] - table['test']).round(3)
print(table.to_string(index=False))
print()
print('An unlimited tree reaches a perfect training score by giving')
print('almost every row its own leaf, and the test score falls.')
Open this example in the editor and change it. Nothing you type is saved.
Random forest
Fit many trees, each on a bootstrap sample of the rows and a random subset of the features at each split, and average their votes. The individual trees are overfitted; their errors are partly independent, so averaging cancels them.
import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, roc_auc_score
rng = np.random.default_rng(7)
n = 3000
X = rng.normal(size=(n, 8))
y = ((X[:, 0] > 0) == (X[:, 1] > 0)).astype(int) # an interaction
y = np.where(rng.random(n) < 0.12, 1 - y, y) # plus noise
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)
single = DecisionTreeClassifier(random_state=0).fit(X_tr, y_tr)
print(f'one unlimited tree: test {accuracy_score(y_te, single.predict(X_te)):.3f} '
f'train {accuracy_score(y_tr, single.predict(X_tr)):.3f}')
rows = []
for trees in [1, 5, 25, 100, 400]:
forest = RandomForestClassifier(n_estimators=trees, random_state=0).fit(X_tr, y_tr)
rows.append({'trees': trees,
'test_acc': round(accuracy_score(y_te, forest.predict(X_te)), 3),
'auc': round(roc_auc_score(y_te, forest.predict_proba(X_te)[:, 1]), 3)})
print()
print(pd.DataFrame(rows).to_string(index=False))
print()
print('More trees never hurts accuracy - it only costs time. That is')
print('why n_estimators is not really a tuning parameter: set it as')
print('high as you can afford.')
Open this example in the editor and change it. Nothing you type is saved.
The hyperparameters that matter, in order:
n_estimators- more is better, with diminishing returns after a few hundred. Not a risk of overfitting.max_depth/min_samples_leaf- the main controls on overfitting.min_samples_leaf=5or more is a good default on noisy data.max_features- how many features each split may consider. Lower means more decorrelated trees.'sqrt'is the default for classification and is usually right.class_weight='balanced'- for a skewed label.
Gradient boosting
Fit a small tree, look at what it got wrong, fit another tree to the remaining error, and repeat. Each tree is weak; the sum is strong.
import numpy as np
import pandas as pd
import time
from sklearn.ensemble import (RandomForestClassifier, GradientBoostingClassifier,
HistGradientBoostingClassifier)
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
rng = np.random.default_rng(11)
n = 6000
X = rng.normal(size=(n, 10))
signal = (1.1 * X[:, 0] - 0.8 * X[:, 1]
+ 1.4 * (X[:, 2] > 0) * (X[:, 3] > 0) # interaction
+ 0.9 * np.abs(X[:, 4])) # non-linearity
y = (signal + rng.normal(0, 1.0, n) > 1.0).astype(int)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)
models = {
'logistic regression': Pipeline([('s', StandardScaler()),
('m', LogisticRegression(max_iter=1000))]),
'random forest': RandomForestClassifier(n_estimators=300, min_samples_leaf=3,
random_state=0),
'gradient boosting': GradientBoostingClassifier(random_state=0),
'hist gradient boosting': HistGradientBoostingClassifier(random_state=0),
}
for name, model in models.items():
start = time.perf_counter()
model.fit(X_tr, y_tr)
seconds = time.perf_counter() - start
auc = roc_auc_score(y_te, model.predict_proba(X_te)[:, 1])
print(f'{name:<24} AUC {auc:.3f} {seconds:6.2f}s')
print()
print('With interactions and non-linearity present, the boosted models')
print('win - which is the usual result on real tabular data.')
Open this example in the editor and change it. Nothing you type is saved.
HistGradientBoostingClassifier is the one to reach for in practice: it bins the features, which makes it far faster on large data, and it handles missing values natively with no imputer.
The three hyperparameters that matter, and they trade off:
import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score
rng = np.random.default_rng(13)
n = 5000
X = rng.normal(size=(n, 8))
y = ((1.0 * X[:, 0] + 1.2 * (X[:, 1] > 0) * (X[:, 2] > 0)
+ rng.normal(0, 1.0, n)) > 0.7).astype(int)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)
rows = []
for rate, iterations in [(0.3, 50), (0.1, 150), (0.05, 400), (0.5, 400)]:
model = HistGradientBoostingClassifier(
learning_rate=rate, max_iter=iterations, max_leaf_nodes=15,
early_stopping=False, random_state=0).fit(X_tr, y_tr)
rows.append({
'learning_rate': rate, 'max_iter': iterations,
'train_auc': round(roc_auc_score(y_tr, model.predict_proba(X_tr)[:, 1]), 3),
'test_auc': round(roc_auc_score(y_te, model.predict_proba(X_te)[:, 1]), 3),
})
print(pd.DataFrame(rows).to_string(index=False))
print()
print('A low learning rate with many iterations generalises best; a')
print('high rate with many iterations overfits. The two must move in')
print('opposite directions.')
Open this example in the editor and change it. Nothing you type is saved.
Use early stopping rather than guessing the iteration count:
import numpy as np
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score
rng = np.random.default_rng(17)
n = 5000
X = rng.normal(size=(n, 8))
y = ((X[:, 0] + 1.2 * (X[:, 1] > 0) * (X[:, 2] > 0) + rng.normal(0, 1, n)) > 0.7).astype(int)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)
model = HistGradientBoostingClassifier(
learning_rate=0.05, max_iter=2000, early_stopping=True,
validation_fraction=0.15, n_iter_no_change=25, random_state=0).fit(X_tr, y_tr)
print(f'max_iter asked for: 2000')
print(f'iterations used: {model.n_iter_}')
print(f'test AUC: {roc_auc_score(y_te, model.predict_proba(X_te)[:, 1]):.3f}')
Open this example in the editor and change it. Nothing you type is saved.
Feature importance, carefully
import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestClassifier
from sklearn.inspection import permutation_importance
from sklearn.model_selection import train_test_split
rng = np.random.default_rng(19)
n = 3000
# Three features matter. Two do not: one is pure noise, and one is a
# high-cardinality id, which impurity importance loves.
useful_a = rng.normal(0, 1, n)
useful_b = rng.integers(0, 3, n)
useful_c = rng.normal(0, 1, n)
noise = rng.normal(0, 1, n)
row_id = np.arange(n)
y = ((1.3 * useful_a + 0.8 * useful_b + 0.6 * useful_c
+ rng.normal(0, 1.0, n)) > 1.0).astype(int)
X = pd.DataFrame({'useful_a': useful_a, 'useful_b': useful_b,
'useful_c': useful_c, 'noise': noise, 'row_id': row_id})
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)
forest = RandomForestClassifier(n_estimators=300, random_state=0).fit(X_tr, y_tr)
impurity = pd.Series(forest.feature_importances_, index=X.columns)
permutation = permutation_importance(forest, X_te, y_te, n_repeats=10, random_state=0)
table = pd.DataFrame({
'impurity': impurity.round(3),
'permutation': pd.Series(permutation.importances_mean, index=X.columns).round(4),
}).sort_values('impurity', ascending=False)
print(table.to_string())
print()
print('Impurity importance ranks row_id and noise far too high: a')
print('continuous feature offers more places to split, so it wins')
print('ties by accident. Permutation importance, computed on held-out')
print('data, puts them near zero where they belong.')
Open this example in the editor and change it. Nothing you type is saved.
Prefer permutation importance, and compute it on the test set. Impurity importance is biased towards high-cardinality features and is computed on training data, which means a feature the model memorised looks important.
Both share a limitation: with two correlated features, importance is split between them arbitrarily, and dropping either one individually costs nothing.
Choosing
import pandas as pd
print(pd.DataFrame([
{'situation': 'tabular data, mixed types, interactions likely',
'reach for': 'HistGradientBoosting, then a random forest'},
{'situation': 'the relationship looks linear',
'reach for': 'logistic or linear regression'},
{'situation': 'the decision must be explainable line by line',
'reach for': 'one shallow tree, or a linear model'},
{'situation': 'inputs will go outside the training range',
'reach for': 'a linear model - trees cannot extrapolate'},
{'situation': 'images, audio, text',
'reach for': 'a neural network; trees are the wrong family'},
{'situation': 'a few hundred rows',
'reach for': 'a linear model with regularisation'},
]).to_string(index=False))
Open this example in the editor and change it. Nothing you type is saved.
The extrapolation limit is worth seeing, because it is absolute:
import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
rng = np.random.default_rng(23)
X = rng.uniform(0, 10, (600, 1))
y = 3 * X[:, 0] + rng.normal(0, 1, 600)
forest = RandomForestRegressor(n_estimators=200, random_state=0).fit(X, y)
linear = LinearRegression().fit(X, y)
print(f'{"input":>8} {"truth":>8} {"forest":>8} {"linear":>8}')
for value in [5.0, 9.9, 15.0, 40.0]:
print(f'{value:>8.1f} {3 * value:>8.1f} '
f'{forest.predict([[value]])[0]:>8.1f} {linear.predict([[value]])[0]:>8.1f}')
print()
print('The forest flat-lines at the edge of what it saw. It cannot do')
print('otherwise: a leaf holds an average of training rows.')
Open this example in the editor and change it. Nothing you type is saved.
When it goes wrong
| Symptom | Cause |
|---|---|
| Perfect training score, poor test | An unlimited tree; set max_depth or min_samples_leaf |
| Predictions flat-line outside a range | Trees cannot extrapolate; use a linear model |
| An id column is the top feature | Impurity importance bias; use permutation |
| Two correlated features both look useless | Importance is split between them |
| More trees made no difference | Expected; n_estimators saturates |
| Boosting overfitted badly | Learning rate too high for the iteration count |
| Training takes minutes | Use HistGradientBoosting instead |
| The minority class is never predicted | Set class_weight='balanced' or move the threshold |
A check you can run
Fit a random forest and a HistGradientBoostingClassifier on your data with default settings, and compare both with the logistic regression from the previous lesson.
Three numbers, ten minutes. If the trees win by a lot, the relationship has interactions or non-linearity and that is worth knowing. If they do not, ship the linear model - it is faster, explainable, and it extrapolates.