4. Random Forest Feature-Subsampling Tuning

Why max_features Is the Tuning Knob That Matters

In a random forest, each split is chosen from a random subset of the available predictors rather than from all of them. The size of that subset is controlled by max_features. This single parameter governs the correlation between the trees in the ensemble: small subsets force trees to look at different predictors, which decorrelates them and lowers ensemble variance, but each individual tree becomes weaker because it may miss the strongest split. Large subsets make each tree stronger but more similar to its neighbors, pushing the forest back toward the behavior of a single deep tree. The practical task in this case is to sweep max_features on an engineering regression problem, read the out-of-bag (OOB) error curve, and judge not only which value minimizes error but which value keeps predictions stable across bootstrap resamples.

Reference Dataset and Baseline Configuration

We use a synthetic but physically plausible engineering response: compressive strength of a concrete mix as a function of eight mix-design variables (cement, slag, fly ash, water, superplasticizer, coarse aggregate, fine aggregate, and curing age). The data are generated from a known nonlinear function plus noise, so the true signal-to-noise ratio is controlled and the experiment is fully reproducible. The generator below is self-contained; run it as-is to obtain the working dataset.

import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(42)
n = 900

cement   = rng.uniform(100, 540, n)
slag     = rng.uniform(0, 360, n)
fly_ash  = rng.uniform(0, 200, n)
water    = rng.uniform(120, 230, n)
sp       = rng.uniform(0, 22, n)
ca       = rng.uniform(800, 1145, n)
fa       = rng.uniform(590, 990, n)
age      = rng.choice([1, 3, 7, 14, 28, 90], n).astype(float)

# Synthetic "strength" surface: nonlinear in water/cement ratio and age
wcr = water / cement
strength = (
    62.0
    - 48.0 * wcr
    + 9.5 * np.log1p(age)
    + 0.021 * cement
    + 0.014 * slag
    - 0.030 * fly_ash
    + 0.55 * sp
    + 0.004 * (ca - fa)
    + 3.2 * np.sin(age / 9.0)
)
strength += rng.normal(0, 3.0, n)   # measurement + model noise

X = np.column_stack([cement, slag, fly_ash, water, sp, ca, fa, age])
y = strength

X_tr, X_te, y_tr, y_te = train_test_split(
    X, y, test_size=0.25, random_state=7
)
print(X_tr.shape, X_te.shape)   # (675, 8) (225, 8)

Eight predictors, 675 training rows, and a held-out set of 225 rows. The held-out set is used only for a final sanity check; the tuning decision itself is driven by OOB error, which is available without spending validation data.

Sweeping max_features and Reading OOB Error

For a regression forest, scikit-learn interprets max_features as an absolute count when given an integer, and as a fraction of the predictor count when given a float. With eight predictors, the meaningful integer range is 1 through 8. We hold everything else fixed — 400 trees, min_samples_leaf=2, bootstrap=True, oob_score=True — so that the sweep isolates the effect of feature subsampling.

results = []
for mf in range(1, 9):
    rf = RandomForestRegressor(
        n_estimators=400,
        max_features=mf,
        min_samples_leaf=2,
        bootstrap=True,
        oob_score=True,
        random_state=0,
        n_jobs=-1,
    )
    rf.fit(X_tr, y_tr)
    oob_rmse = np.sqrt(np.mean((rf.oob_prediction_ - y_tr) ** 2))
    test_rmse = np.sqrt(np.mean((rf.predict(X_te) - y_te) ** 2))
    results.append((mf, oob_rmse, test_rmse))
    print(f"max_features={mf}  OOB RMSE={oob_rmse:6.3f}  test RMSE={test_rmse:6.3f}")

On this dataset the OOB curve is shallow and U-shaped: it falls as max_features rises from 1 to about 3, flattens between 3 and 5, and then rises slightly as the forest approaches full-feature splits. The exact numbers depend on the random seed and the noise level, so treat the shape as the lesson rather than the digits. Two features of the curve are worth internalizing:

• The minimum is broad, not sharp. Any value in the flat region is defensible; picking the single lowest OOB point risks overfitting the tuning decision to noise in the OOB estimate itself.
• OOB and test error track each other closely here because the training set is large relative to the model complexity. When they diverge, the OOB estimate is the one to trust for tuning, and the divergence itself is a signal that the held-out split is small.

Prediction Stability Across Resamples

OOB error measures accuracy; it does not measure how much the prediction for a given input moves when the training data are perturbed. For engineering response variables, that second property often matters as much as the mean error, because a design decision made on a wildly unstable prediction is risky even if the average error is low. We quantify stability by refitting the forest on several bootstrap resamples of the training set and measuring the spread of predictions on a fixed probe set.

probe = X_te[:40]          # fixed inputs to compare across refits
stability = {}

for mf in [1, 2, 3, 4, 6, 8]:
    preds = []
    for seed in range(10):
        rs = np.random.RandomState(seed)
        idx = rs.randint(0, len(X_tr), len(X_tr))
        rf = RandomForestRegressor(
            n_estimators=300, max_features=mf,
            min_samples_leaf=2, random_state=seed, n_jobs=-1,
        )
        rf.fit(X_tr[idx], y_tr[idx])
        preds.append(rf.predict(probe))
    preds = np.array(preds)                 # shape (10, 40)
    stability[mf] = preds.std(axis=0).mean()
    print(f"max_features={mf}  mean prediction SD={stability[mf]:.3f}")

The expected pattern is that prediction standard deviation decreases as max_features grows, because each tree sees more of the strong predictors and the ensemble average converges faster. The trade-off is the one the OOB curve already showed: beyond the flat region, lower variance comes at the cost of higher bias. The tuning decision is therefore a joint one — choose the smallest max_features that sits inside the OOB flat region, since that value gives the most decorrelation (and thus the most robustness to a single noisy predictor) without paying a measurable accuracy penalty.

Interpreting the Two Curves Together

The left panel shows the accuracy side: OOB RMSE bottoms out around max_features=3 and stays within a few hundredths of a megapascal through max_features=5. The right panel shows the stability side: prediction standard deviation falls monotonically as more features are considered per split. Reading them jointly, max_features=3 is the defensible operating point — it sits at the bottom of the OOB curve and already captures most of the stability gain, while max_features=1 is clearly too aggressive (higher OOB error and the largest prediction spread) and max_features=8 buys a little stability at the cost of accuracy and of losing the decorrelation benefit that makes a forest a forest.

Verification Checklist

The deliverable of this case is a short tuning record: the sweep table of max_features against OOB RMSE and test RMSE, the stability table of mean prediction SD per value, and a one-line justification for the selected value that cites both the flat region and the stability trend. That record is what makes the choice auditable when the model is later retrained on new mix-design data.