Model Building¶
Training a model without a baseline is guesswork. The baseline anchors your expectations — it tells you how much of the predictive problem is solved by something trivially simple. Every subsequent model should beat it by a meaningful margin to justify its complexity.
Setup — Rebuild Features from Previous Steps¶
This file continues from Feature Engineering. Run the full pipeline first.
from sklearn.datasets import fetch_california_housing
import pandas as pd
import numpy as np
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
housing = fetch_california_housing(as_frame=True)
df = housing.frame.copy()
# --- Outlier capping ---
def cap_at_percentile(series, upper_pct=99):
upper = series.quantile(upper_pct / 100)
return series.clip(upper=upper)
for col in ['AveRooms', 'AveBedrms', 'AveOccup']:
df[col] = cap_at_percentile(df[col], upper_pct=99)
# --- Interaction features ---
df['rooms_per_person'] = df['AveRooms'] / df['AveOccup']
df['bedrooms_ratio'] = df['AveBedrms'] / df['AveRooms']
df['income_per_room'] = df['MedInc'] / df['AveRooms']
# --- Log transforms ---
df['log_population'] = np.log1p(df['Population'])
df['log_ave_occup'] = np.log1p(df['AveOccup'])
# --- Geographic features ---
def simplified_coast_distance(lat, lon):
coast_points = [
(37.8, -122.4), (36.6, -121.9), (34.4, -119.7),
(33.9, -118.4), (32.7, -117.2),
]
return min(np.sqrt((lat - c_lat)**2 + (lon - c_lon)**2)
for c_lat, c_lon in coast_points)
df['coast_distance'] = df.apply(
lambda row: simplified_coast_distance(row['Latitude'], row['Longitude']), axis=1
)
kmeans = KMeans(n_clusters=6, random_state=42, n_init=10)
df['geo_cluster'] = kmeans.fit_predict(df[['Latitude', 'Longitude']].values)
# --- Feature matrix and target ---
feature_cols = [
'MedInc', 'HouseAge', 'AveRooms', 'AveBedrms', 'Population',
'AveOccup', 'Latitude', 'Longitude',
'rooms_per_person', 'bedrooms_ratio', 'income_per_room',
'log_population', 'log_ave_occup',
'coast_distance', 'geo_cluster'
]
X = df[feature_cols].copy()
y = df['MedHouseVal'].copy()
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42
)
# Scaled versions for linear models
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
print(f"Training set: {X_train.shape}") # (16512, 15)
print(f"Test set: {X_test.shape}") # (4128, 15)
Helper — Evaluation Function¶
Define this once and call it for every model.
def evaluate_model(model_name, y_train, y_train_pred, y_test, y_test_pred):
"""Print MAE, RMSE and R² for train and test sets."""
train_mae = mean_absolute_error(y_train, y_train_pred)
test_mae = mean_absolute_error(y_test, y_test_pred)
train_rmse = np.sqrt(mean_squared_error(y_train, y_train_pred))
test_rmse = np.sqrt(mean_squared_error(y_test, y_test_pred))
train_r2 = r2_score(y_train, y_train_pred)
test_r2 = r2_score(y_test, y_test_pred)
print(f"\n{'='*50}")
print(f" {model_name}")
print(f"{'='*50}")
print(f" {'Metric':<12} {'Train':>10} {'Test':>10}")
print(f" {'-'*32}")
print(f" {'MAE':<12} {train_mae:>10.4f} {test_mae:>10.4f}")
print(f" {'RMSE':<12} {train_rmse:>10.4f} {test_rmse:>10.4f}")
print(f" {'R²':<12} {train_r2:>10.4f} {test_r2:>10.4f}")
return {
'model': model_name,
'train_mae': train_mae, 'test_mae': test_mae,
'train_rmse': train_rmse, 'test_rmse': test_rmse,
'train_r2': train_r2, 'test_r2': test_r2,
}
results = [] # collect all model results for the comparison table
Baseline 1 — Predict the Training Mean¶
The simplest possible model: predict the same value for every house. This tells you the floor — what you get from zero predictive effort.
train_mean = y_train.mean()
mean_preds_train = np.full(len(y_train), train_mean)
mean_preds_test = np.full(len(y_test), train_mean)
results.append(evaluate_model(
'Baseline (Mean)',
y_train, mean_preds_train,
y_test, mean_preds_test
))
# ==================================================
# Baseline (Mean)
# ==================================================
# Metric Train Test
# --------------------------------
# MAE 0.8617 0.8614
# RMSE 1.1543 1.1555
# R² 0.0000 0.0000
An RMSE of 1.16 (in units of $100k) means that predicting the average gives you an error of about $116,000 on average. R² = 0.0 by definition — the mean model explains none of the variance.
Baseline 2 — Predict the Training Median¶
The median is more robust than the mean when the target is skewed or capped.
train_median = y_train.median()
median_preds_train = np.full(len(y_train), train_median)
median_preds_test = np.full(len(y_test), train_median)
results.append(evaluate_model(
'Baseline (Median)',
y_train, median_preds_train,
y_test, median_preds_test
))
# ==================================================
# Baseline (Median)
# ==================================================
# MAE 0.7930 0.7927
# RMSE 1.1580 1.1592
# R² -0.0063 -0.0065
MAE drops slightly versus the mean baseline because MAE is minimized by the median, but RMSE is similar. The negative R² is expected — neither baseline model explains any variance.
Info
These two baselines answer a specific question: how hard is this problem? An RMSE of 1.16 ($116k) for a naive guess tells you the target has meaningful variance. A model with R² = 0.80 has reduced that to about 0.52 RMSE ($52k). Without a baseline, you have no way to know if 0.52 is impressive or disappointing.
Linear Regression¶
Linear regression on scaled features. This is the simplest model that actually learns from the data.
from sklearn.linear_model import LinearRegression
lr = LinearRegression()
lr.fit(X_train_scaled, y_train)
lr_train_pred = lr.predict(X_train_scaled)
lr_test_pred = lr.predict(X_test_scaled)
results.append(evaluate_model(
'Linear Regression',
y_train, lr_train_pred,
y_test, lr_test_pred
))
# ==================================================
# Linear Regression
# ==================================================
# MAE 0.5339 0.5354
# RMSE 0.7285 0.7326
# R² 0.6019 0.5993
# Inspect the most influential coefficients
coef_df = pd.DataFrame({
'feature': feature_cols,
'coefficient': lr.coef_
}).sort_values('coefficient', key=abs, ascending=False)
print(coef_df.head(8).to_string(index=False))
# feature coefficient
# MedInc 0.4512
# Longitude -0.3601
# Latitude -0.3148
# rooms_per_person 0.1823
# income_per_room 0.1347
# HouseAge 0.1021
# bedrooms_ratio -0.0934
# log_ave_occup -0.0879
Info
Coefficients are on the standardized scale, so they are comparable. MedInc has the largest positive coefficient — a one-standard-deviation increase in income is associated with a 0.45 unit ($45k) increase in predicted house value. The large negative coefficients on Latitude and Longitude reflect California's north-south and inland-coastal gradients, but a linear model cannot capture their interaction — it can only add them independently.
Ridge Regression¶
Ridge adds L2 regularization, which penalizes large coefficients. This helps when features are correlated (as MedInc, income_per_room, and rooms_per_person are).
from sklearn.linear_model import Ridge
ridge = Ridge(alpha=1.0)
ridge.fit(X_train_scaled, y_train)
ridge_train_pred = ridge.predict(X_train_scaled)
ridge_test_pred = ridge.predict(X_test_scaled)
results.append(evaluate_model(
'Ridge (alpha=1.0)',
y_train, ridge_train_pred,
y_test, ridge_test_pred
))
# ==================================================
# Ridge (alpha=1.0)
# ==================================================
# MAE 0.5337 0.5351
# RMSE 0.7283 0.7322
# R² 0.6021 0.5997
The improvement over plain linear regression is small but consistent. Ridge has slightly lower test RMSE because it shrinks the most volatile coefficients, reducing variance at the cost of a tiny increase in bias.
Tip
The small gap between Linear Regression and Ridge here suggests that multicollinearity is not a severe problem in this dataset. The bigger gains come from non-linear models, not from regularization. Ridge is worth trying, but do not expect dramatic improvement on data that is already reasonably well-conditioned.
Random Forest Regressor¶
Random Forest builds many decision trees on bootstrapped samples and averages their predictions. It handles non-linear relationships and feature interactions naturally.
from sklearn.ensemble import RandomForestRegressor
rf = RandomForestRegressor(n_estimators=100, random_state=42, n_jobs=-1)
rf.fit(X_train, y_train) # tree models use unscaled features
rf_train_pred = rf.predict(X_train)
rf_test_pred = rf.predict(X_test)
results.append(evaluate_model(
'Random Forest (n=100)',
y_train, rf_train_pred,
y_test, rf_test_pred
))
# ==================================================
# Random Forest (n=100)
# ==================================================
# MAE 0.1734 0.3290
# RMSE 0.2612 0.5115
# R² 0.9489 0.8036
The jump from R² ≈ 0.60 (linear) to R² ≈ 0.80 (Random Forest) on the test set is the payoff for switching to a non-linear model. Notice, however, that training RMSE (0.26) is much lower than test RMSE (0.51). That gap is overfitting — the forest has memorized training samples.
Warning
A train R² of 0.95 and a test R² of 0.80 is a sign of moderate overfitting. The model is learning some patterns that are specific to training instances rather than generalizable structure. This is common with default Random Forest settings. You will address this with hyperparameter tuning below.
Gradient Boosting Regressor¶
Gradient Boosting builds trees sequentially, where each tree corrects the residual errors of the previous one. It tends to be more accurate than Random Forest but slower to train.
from sklearn.ensemble import GradientBoostingRegressor
gbr = GradientBoostingRegressor(
n_estimators=300,
learning_rate=0.05,
max_depth=5,
subsample=0.8,
random_state=42
)
gbr.fit(X_train, y_train)
gbr_train_pred = gbr.predict(X_train)
gbr_test_pred = gbr.predict(X_test)
results.append(evaluate_model(
'Gradient Boosting',
y_train, gbr_train_pred,
y_test, gbr_test_pred
))
# ==================================================
# Gradient Boosting
# ==================================================
# MAE 0.2761 0.3143
# RMSE 0.3827 0.4549
# R² 0.8899 0.8453
Gradient Boosting achieves the best test RMSE so far (0.455) with a tighter train-test gap than Random Forest, meaning it generalizes better despite fitting the training data less aggressively.
Model Comparison Table¶
import matplotlib.pyplot as plt
comparison = pd.DataFrame(results)
print(comparison[['model','train_mae','test_mae','train_rmse','test_rmse','train_r2','test_r2']].round(4).to_string(index=False))
# model train_mae test_mae train_rmse test_rmse train_r2 test_r2
# Baseline (Mean) 0.8617 0.8614 1.1543 1.1555 0.0000 0.0000
# Baseline (Median) 0.7930 0.7927 1.1580 1.1592 -0.0063 -0.0065
# Linear Regression 0.5339 0.5354 0.7285 0.7326 0.6019 0.5993
# Ridge (alpha=1.0) 0.5337 0.5351 0.7283 0.7322 0.6021 0.5997
# Random Forest (n=100) 0.1734 0.3290 0.2612 0.5115 0.9489 0.8036
# Gradient Boosting 0.2761 0.3143 0.3827 0.4549 0.8899 0.8453
Success
Gradient Boosting is the best single model: lowest test MAE (0.314, roughly $31k average error) and highest test R² (0.845). Its train-test gap is also smaller than Random Forest's, indicating better generalization.
Hyperparameter Tuning — Random Forest with RandomizedSearchCV¶
RandomizedSearchCV samples random combinations from the parameter grid rather than exhaustively trying every combination. This is much faster than GridSearchCV while still finding a near-optimal configuration.
from sklearn.model_selection import RandomizedSearchCV
param_dist = {
'n_estimators': [100, 200, 300],
'max_depth': [None, 10, 20],
'min_samples_leaf': [1, 2, 5],
'max_features': ['sqrt', 0.5],
}
rf_base = RandomForestRegressor(random_state=42, n_jobs=-1)
rf_search = RandomizedSearchCV(
rf_base,
param_distributions=param_dist,
n_iter=20,
cv=3,
scoring='neg_root_mean_squared_error',
random_state=42,
n_jobs=-1,
verbose=1
)
rf_search.fit(X_train, y_train)
print("Best parameters:", rf_search.best_params_)
# Best parameters: {'n_estimators': 300, 'min_samples_leaf': 2,
# 'max_features': 'sqrt', 'max_depth': 20}
print(f"Best CV RMSE: {-rf_search.best_score_:.4f}")
# Best CV RMSE: 0.4921
best_rf = rf_search.best_estimator_
best_rf_train_pred = best_rf.predict(X_train)
best_rf_test_pred = best_rf.predict(X_test)
results.append(evaluate_model(
'Random Forest (Tuned)',
y_train, best_rf_train_pred,
y_test, best_rf_test_pred
))
# ==================================================
# Random Forest (Tuned)
# ==================================================
# MAE 0.1991 0.3162
# RMSE 0.3107 0.4935
# R² 0.9272 0.8181
The tuned Random Forest closes much of the train-test gap and now competes with Gradient Boosting on test RMSE.
Tip
min_samples_leaf=2 is doing meaningful work here. It prevents leaves from fitting single training instances, which is the main source of overfitting in the default Random Forest. When you see a large train-test RMSE gap, try increasing min_samples_leaf before adding more trees.
The Ceiling Effect — Predicted vs. Actual¶
One pattern that no amount of tuning will fully fix is the ceiling at 5.0.
import matplotlib.pyplot as plt
plt.figure(figsize=(7, 6))
plt.scatter(y_test, best_rf_test_pred, alpha=0.15, s=5, color='steelblue')
plt.axline((0, 0), slope=1, color='red', linewidth=1.5, label='Perfect prediction')
plt.axvline(5.0, color='orange', linestyle='--', linewidth=1.2, label='Ceiling at 5.0')
plt.xlabel('Actual MedHouseVal ($100k)')
plt.ylabel('Predicted MedHouseVal ($100k)')
plt.title('Predicted vs. Actual — Tuned Random Forest')
plt.legend()
plt.tight_layout()
plt.show()
The orange dashed line reveals the problem. At actual = 5.0 there is a horizontal cluster of points: hundreds of blocks where the true value was capped at $500k but the actual market value was higher. The model predicts below 5.0 for many of these — not because it is wrong, but because the training data never showed it what happens above the ceiling. No model trained on this data can accurately predict those blocks.
Warning
This is a data quality issue, not a modeling issue. The ceiling at $500k was imposed by the 1990 census data collection process. Any model you train on this data inherits this blind spot. Document this limitation clearly before deploying any predictions to stakeholders. The model is unreliable for high-value properties — precisely the ones stakeholders often care most about.