Random Forest: Forest Cover Type Project

This is the Random Forest capstone: a full pipeline from raw data to deployment-ready model on Forest Cover Type — 580k samples, 54 features, 7 classes. It demonstrates RF's scalability, its calibration behavior, and the realistic gap between a tuned RF and a single decision tree.

Anchor dataset: Forest Cover Type (sklearn's fetch_covtype) — cartographic features to predict which of 7 tree species covers a 30×30m forest patch.

from sklearn.datasets import fetch_covtype
from sklearn.model_selection import train_test_split
import pandas as pd
import numpy as np

data = fetch_covtype()
X_full, y_full = data.data, data.target - 1  # convert to 0-indexed classes

# 10% subset for speed (58,101 samples)
X, _, y, _ = train_test_split(X_full, y_full, train_size=0.10,
                               random_state=42, stratify=y_full)

Step 1: Dataset Overview

print(f"Shape: {X.shape}")
print(f"Continuous features: {X[:, :10].shape[1]} (elevation, slope, distances, hillshade)")
print(f"Binary features:     {X[:, 10:].shape[1]} (wilderness area + soil type OHE)")
print(f"\nClass distribution:")
for cls in range(7):
    count = (y == cls).sum()
    print(f"  Class {cls} ({data.target_names[cls]}): {count:5d} ({count/len(y)*100:.1f}%)")

54 features = 10 continuous (elevation, slope, horizontal/vertical distance to water, fire points, roadways, 3 hillshade values) + 44 binary (4 wilderness area + 40 soil type one-hot encoded). Heavy class imbalance: Classes 1 and 0 make up 84% of samples.

Step 2: Train/Test Split

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y
)
print(f"Train: {X_train.shape}, Test: {X_test.shape}")

Step 3: Baseline — Single Decision Tree

from sklearn.tree import DecisionTreeClassifier

dt = DecisionTreeClassifier(random_state=42)
dt.fit(X_train, y_train)
print(f"Decision Tree: Train={dt.score(X_train, y_train):.4f}, Test={dt.score(X_test, y_test):.4f}")
print(f"Depth: {dt.get_depth()}, Leaves: {dt.get_n_leaves()}")

The tree memorizes all 46,480 training samples (nearly one leaf per sample at depth 43). Test accuracy of 84% shows the expected overfitting gap.

Step 4: Random Forest — Default Settings

from sklearn.ensemble import RandomForestClassifier

rf_default = RandomForestClassifier(n_estimators=100, oob_score=True,
                                     random_state=42, n_jobs=-1)
rf_default.fit(X_train, y_train)
print(f"RF (default): Train={rf_default.score(X_train, y_train):.4f}, Test={rf_default.score(X_test, y_test):.4f}")
print(f"OOB accuracy: {rf_default.oob_score_:.4f}")

100 trees with default settings: test accuracy jumps from 0.840 (single tree) to 0.928 — an 88-point improvement. OOB (0.925) closely tracks test (0.928), confirming reliable free validation.

Step 5: Hyperparameter Tuning

GridSearchCV over 54 features × 46k samples is expensive. RandomizedSearchCV samples from the space instead:

from sklearn.model_selection import RandomizedSearchCV

param_dist = {
    'n_estimators':    [50, 100, 200],
    'max_depth':       [None, 10, 20, 30],
    'max_features':    ['sqrt', 'log2', 0.5],
    'min_samples_leaf':  [1, 2, 5],
    'min_samples_split': [2, 5, 10],
}
rs = RandomizedSearchCV(
    RandomForestClassifier(random_state=42, n_jobs=-1, oob_score=True),
    param_dist, n_iter=20, cv=5, scoring='accuracy',
    random_state=42, n_jobs=-1
)
rs.fit(X_train, y_train)
print(f"Best params:    {rs.best_params_}")
print(f"Best CV acc:    {rs.best_score_:.4f}")
print(f"Test accuracy:  {rs.best_estimator_.score(X_test, y_test):.4f}")

Tuned RF (n=200, min_samples_leaf=2) adds +0.0076 over the default — a small but real gain. min_samples_leaf=2 prevents leaves with single samples, which reduces overfitting without restricting tree depth.

Step 6: Full Evaluation — Classification Report

from sklearn.metrics import classification_report, confusion_matrix

best_rf = rs.best_estimator_
y_pred = best_rf.predict(X_test)
class_names = [f"Class_{i}" for i in range(7)]

print(classification_report(y_test, y_pred, target_names=class_names))

Classes 0 and 1 (dominant at 84% combined) have F1 ≥ 0.95. Rarer classes (Class 4 Aspen: F1=0.84, Class 6 Krummholz: F1=0.89) have slightly lower recall — fewer training examples to learn from.

Step 7: Feature Importance

importances = pd.Series(best_rf.feature_importances_,
                         index=data.feature_names).sort_values(ascending=False)
print("Top 10 features:")
print(importances.head(10).round(4))

Elevation dominates at 23% — makes ecological sense, as different tree species occupy different altitude bands. All top 10 are continuous features; the 44 binary soil-type features share the remaining ~25% importance, each contributing under 1%.

Step 8: Confusion Matrix

cm = confusion_matrix(y_test, y_pred)
print("Confusion Matrix (rows=true, cols=predicted):")
print(cm)

Main error pattern: Class 1 (Lodgepole Pine) occasionally misclassified as Class 0 (Spruce/Fir) and vice versa — 195+148 errors. These species share overlapping elevation ranges and are the dataset's hardest distinction. Class 6 (Krummholz) — 65 test samples — achieves perfect recall (63/65=96.9%) despite being rare.

Step 9: Model Calibration

A well-calibrated model is more confident when it's right than when it's wrong:

y_prob = best_rf.predict_proba(X_test)

correct_mask = (y_pred == y_test)
max_prob_correct = y_prob[correct_mask].max(axis=1)
max_prob_wrong   = y_prob[~correct_mask].max(axis=1)

print(f"Mean confidence when CORRECT: {max_prob_correct.mean():.4f}")
print(f"Mean confidence when WRONG:   {max_prob_wrong.mean():.4f}")
print(f"Fraction >90% confident when correct: {(max_prob_correct > 0.9).mean():.4f}")
print(f"Fraction >90% confident when wrong:   {(max_prob_wrong > 0.9).mean():.4f}")

When the model is correct, average confidence is 89%. When it's wrong, it drops to 56% — barely above the majority-class rate. The RF's probability estimates (computed as fraction of trees predicting each class) are practically calibrated: high confidence predictions are reliable.

Step 10: Save and Load Model

import joblib

joblib.dump(best_rf, 'forest_cover_rf.pkl')

# Verify round-trip
loaded = joblib.load('forest_cover_rf.pkl')
sample_pred = loaded.predict(X_test[:5])
print(f"Sample predictions: {sample_pred}")
print(f"True labels:        {y_test[:5].tolist()}")

All 5 match. joblib is faster than pickle for sklearn models with large numpy arrays (trees store split thresholds as float arrays).

Performance Summary

Test Your Understanding

1. The default RF achieves Train=0.999 (not 1.0). Why doesn't a Random Forest of decision trees with no max_depth reach 100% training accuracy? Each tree is trained on a bootstrap sample — what fraction of training samples does each tree never see, and how does this affect training accuracy?

2. RandomizedSearchCV with n_iter=20 explores 20 of the 3×4×3×3×3=324 possible parameter combinations. If the true best combination was one of the 304 not tried, how would you know? What evidence could you use — other than running more iterations — to gain confidence that the found parameters are near-optimal?

3. Classes 0 and 1 (Spruce/Fir and Lodgepole Pine) share 195+148=343 errors between each other. These two species dominate the dataset (84% combined). Would SMOTE-oversampling the rare classes (2–6) improve or worsen the confusion between classes 0 and 1? Explain why.

4. The calibration check shows confidence 0.89 when correct vs 0.56 when wrong. Random Forest's probabilities are computed as the fraction of trees voting for each class. If 90 out of 100 trees vote for class 1: probability = 0.90. If you wanted to use RF probabilities in a downstream cost-sensitive system (where misclassifying class 5 as class 1 costs 10× more), what modification to the decision rule would you make?

5. Elevation accounts for 23% of feature importance but the dataset has 44 binary soil-type features. Despite having 44× more binary features than Elevation, their combined importance is ~25%. Why does Elevation dominate MDI so completely, and is this consistent or problematic given what we know about MDI's bias toward high-cardinality features?