Evaluation and Interpretation

A model that performs well in training means nothing until you can explain what "well" means in terms a stakeholder can act on. This section converts raw metrics into business language, visualizes where the model succeeds and fails, and identifies the structural limitations that every user of this model needs to understand.

All code here assumes the tuned Random Forest from Model Building is your primary model. Substitute any other trained model as needed.

Setup — Load Model and Predictions

# Assumes all feature engineering and model training steps have been run.
# best_rf, X_train, X_test, y_train, y_test are available from model-building.md.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

y_pred = best_rf.predict(X_test)
residuals = y_test.values - y_pred

Section 1 — Metrics on the Test Set

mae  = mean_absolute_error(y_test, y_pred)
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
r2   = r2_score(y_test, y_pred)

# MAPE — exclude near-zero targets to avoid division instability
mask = y_test > 0.1
mape = np.mean(np.abs((y_test[mask].values - y_pred[mask]) / y_test[mask].values)) * 100

print(f"MAE:   {mae:.4f}  (${mae * 100_000:,.0f})")
print(f"RMSE:  {rmse:.4f}  (${rmse * 100_000:,.0f})")
print(f"MAPE:  {mape:.2f}%")
print(f"R²:    {r2:.4f}")

# MAE:   0.3162  ($31,620)
# RMSE:  0.4935  ($49,350)
# MAPE:  17.84%
# R²:    0.8181

Interpreting Each Metric in Business Terms

Metric	Value	What it means
MAE = 0.316	$31,600	On average, predictions are off by $31,600. Half of predictions will be closer; half will be further.
RMSE = 0.494	$49,400	The typical error, penalizing large mistakes more heavily. A prediction that is $150k wrong hurts RMSE disproportionately.
MAPE = 17.8%	—	On average, predictions are 17.8% off relative to the true value. A $300k home might be predicted as anywhere from $247k to $354k.
R² = 0.818	—	The model explains 81.8% of the variance in house prices across the test set. The remaining 18.2% is driven by factors not in the dataset.

Info

MAE and RMSE are in the same units as the target (hundreds of thousands of dollars). Divide by the target mean (~$207k) to get a sense of relative scale. An MAE of $31k on a $207k average price is about 15% relative error — reasonable for block-level prediction from census data, but too wide for individual home appraisals.

Tip

Report MAE to business stakeholders — it is the most intuitive. Report RMSE in technical contexts when large errors are especially costly. Both tell a consistent story here.

Section 2 — Residuals Plot

A residuals plot (residuals vs. predicted values) reveals whether errors are random or systematic.

plt.figure(figsize=(9, 5))
plt.scatter(y_pred, residuals, alpha=0.15, s=5, color='steelblue')
plt.axhline(0, color='red', linewidth=1.5, linestyle='--', label='Zero error')
plt.axhline(residuals.mean(), color='orange', linewidth=1, linestyle=':',
            label=f'Mean residual = {residuals.mean():.3f}')
plt.xlabel('Predicted MedHouseVal ($100k)')
plt.ylabel('Residual (Actual − Predicted)')
plt.title('Residuals vs. Predicted Values')
plt.legend()
plt.tight_layout()
plt.show()

print(f"Mean residual:   {residuals.mean():.4f}")
print(f"Std of residuals: {residuals.std():.4f}")

# Mean residual:    0.0193
# Std of residuals: 0.4929

Two patterns should stand out:

1. Funnel shape at the right. As predicted values approach 5.0 ($500k), the residuals become more negative. The model predicts values just below 5.0, but the true values for high-priced blocks are also 5.0 (the ceiling). This creates systematic underestimation — a cluster of negative residuals on the right side of the plot.

2. Roughly zero-centered elsewhere. For predicted values below 4.0, residuals scatter around zero without a strong trend, which is the behavior you want from a well-specified model.

Warning

The systematic underestimation at high predicted values is not random noise — it is a structural bias introduced by the data ceiling. It cannot be corrected through better algorithms. It can only be addressed by using different data (actual transaction prices rather than census-capped medians) or by treating high-value blocks as a separate prediction problem.

Section 3 — Predicted vs. Actual Scatter

plt.figure(figsize=(7, 6))
plt.scatter(y_test, y_pred, alpha=0.15, s=5, color='steelblue', label='Predictions')
plt.axline((0, 0), slope=1, color='red', linewidth=1.5, label='Perfect prediction (y=x)')
plt.axvline(5.0, color='orange', linestyle='--', linewidth=1.2, label='Ceiling at actual = 5.0')
plt.xlabel('Actual MedHouseVal ($100k)')
plt.ylabel('Predicted MedHouseVal ($100k)')
plt.title('Predicted vs. Actual — Tuned Random Forest')
plt.legend(loc='upper left')
plt.xlim(0, 6.5)
plt.ylim(0, 6.5)
plt.tight_layout()
plt.show()

The vertical cluster at actual = 5.0 is the ceiling effect made visible. Hundreds of test blocks have a true value of exactly 5.0, but the model predicts a range of values below 5.0. These blocks were capped by the census collection process — they could have been worth $600k, $700k, or more. The model has no way to know.

For blocks below the ceiling, predictions scatter tightly around the red line, which is the behavior you want.

Success

Below the $500k ceiling, the model is reliable. Above it, it is structurally blind. Communicate this boundary clearly to anyone who uses these predictions for decision-making.

Section 4 — Feature Importance

Random Forest records how much each feature reduces impurity (variance) across all splits in all trees.

importances = best_rf.feature_importances_
feat_importance = pd.Series(importances, index=feature_cols).sort_values(ascending=True)

# Plot top 10
top10 = feat_importance.tail(10)

plt.figure(figsize=(8, 5))
colors = ['#0D9488' if imp > 0.05 else '#64B5C0' for imp in top10]
top10.plot(kind='barh', color=colors)
plt.xlabel('Mean Decrease in Impurity (Importance)')
plt.title('Feature Importance — Tuned Random Forest (Top 10)')
plt.tight_layout()
plt.show()

print(feat_importance.tail(10).round(4).sort_values(ascending=False))

# MedInc            0.4087
# Latitude          0.1423
# Longitude         0.1351
# income_per_room   0.0713
# AveOccup          0.0581
# rooms_per_person  0.0432
# HouseAge          0.0389
# coast_distance    0.0321
# AveRooms          0.0258
# log_population    0.0179

MedInc accounts for roughly 40% of all variance reduction in the forest — confirming the correlation analysis from EDA. Geographic coordinates (Latitude, Longitude) together contribute another 27%, reflecting the strong coastal and regional price gradients in California. The engineered features (income_per_room, rooms_per_person) contribute modestly but consistently, justifying the feature engineering step.

Warning

Feature importance from Random Forest measures how often a feature is used in splits, weighted by the variance reduction it causes. It does not measure causal influence. MedInc being at the top does not mean income causes high prices — both may be caused by a third factor (desirable neighborhood, good schools, proximity to jobs). Use importance for modeling decisions, not causal claims.

Section 5 — Spatial Error Analysis

If prediction errors cluster geographically, it tells you the model is missing a spatial signal — perhaps a neighborhood effect, school district quality, or amenity proximity that is not in the dataset.

# Need lat/lon from the test split — rebuild using the original index
df_test = df.iloc[X_test.index].copy()   # df is the full engineered dataframe
df_test['predicted'] = y_pred
df_test['abs_residual'] = np.abs(y_test.values - y_pred)

plt.figure(figsize=(9, 7))
sc = plt.scatter(
    df_test['Longitude'],
    df_test['Latitude'],
    c=df_test['abs_residual'],
    cmap='YlOrRd',
    alpha=0.6,
    s=4,
    vmin=0,
    vmax=1.5
)
plt.colorbar(sc, label='Absolute Residual ($100k)')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('Spatial Distribution of Prediction Errors')
plt.tight_layout()
plt.show()

# High-error regions
high_error = df_test[df_test['abs_residual'] > 1.0]
print(f"Blocks with error > $100k: {len(high_error)} ({100 * len(high_error) / len(df_test):.1f}%)")
print(high_error[['Latitude', 'Longitude', 'MedHouseVal', 'predicted', 'abs_residual']].head(8).round(3))

Info

High-error blocks concentrate in two areas. First, the coastal zones around Los Angeles, the San Francisco Bay, and the Santa Barbara coast — these are the high-value blocks that hit the $500k ceiling, and the model cannot predict them accurately for the reason discussed above. Second, small clusters in the Central Valley, likely where the block-level census data is noisier due to agricultural or seasonal housing.

Section 6 — Business Interpretation

# What does MAE mean for a typical block?
target_mean = y_test.mean()
print(f"Mean house value in test set: ${target_mean * 100_000:,.0f}")
print(f"MAE as % of mean:             {mae / target_mean * 100:.1f}%")
print(f"R² on test set:               {r2:.3f}")

# Mean house value in test set: $206,800
# MAE as % of mean:             15.3%
# R² on test set:               0.818

For a median-income census block (MedInc ≈ 3.5, roughly $35,000/year household income), the model predicts median home value within approximately $31,600 on average, which represents a 15% relative error. For use cases like regional market monitoring, portfolio valuation, or identifying investment zones, this is sufficient. For use cases like individual home appraisal or mortgage underwriting, a 15% error is too wide.

The model is well-suited for: - Ranking neighborhoods by predicted price tier - Identifying blocks likely to be undervalued relative to their income and geographic characteristics - Estimating aggregate portfolio values across a large number of blocks

The model is not well-suited for: - Pricing individual homes (block-level averages, not individual properties) - Appraising high-value properties (ceiling effect above $500k) - Any application in the current market (data is from 1990)

Section 7 — Limitations

Warning

The $500k ceiling. 4.7% of blocks have a true value of exactly $500k in this dataset. The actual market values of those blocks were unknown at census collection time. The model systematically underestimates these blocks. Any downstream decision that relies on predictions for high-value areas should be treated with extra skepticism.

Warning

1990 vintage data. The California Housing dataset was collected for the 1990 census. California's housing market has changed dramatically since then — values in many coastal areas have increased 5–10x. Relationships learned from 1990 data (income → price, geography → price) may still hold directionally, but the absolute values are not meaningful in the current market without recalibration.

Warning

Block-level aggregates, not individual homes. MedHouseVal is the median value across all homes in a census block. A block may contain a mix of $200k and $800k homes, and the model predicts their median. It tells you nothing about the distribution within a block. Do not use block-level predictions to price individual properties.

Section 8 — What Would You Do Next?

Tip

Add external data sources. School quality ratings (GreatSchools API), walkability scores (Walk Score), distance to transit stops, and crime statistics are all known drivers of residential prices that are absent from this dataset. Adding even one high-quality external signal would likely push R² above 0.90.

Tip

Train a separate model for high-value blocks. Blocks with predicted values above 4.0 behave differently from average blocks. A separate model trained only on high-value blocks — ideally using data without the $500k ceiling — would improve accuracy where the current model is weakest.

Tip

Use quantile regression to model uncertainty. A point prediction of $310k is less useful than an interval: "the true value is likely between $270k and $360k with 80% confidence." Quantile regression (or GradientBoostingRegressor with loss='quantile') produces prediction intervals, which let stakeholders understand model confidence on a per-prediction basis.

Tip

Try XGBoost or LightGBM. These gradient boosting libraries are typically faster and more accurate than scikit-learn's GradientBoostingRegressor. LightGBM in particular handles large datasets efficiently and tends to produce tighter train-test gaps through its leaf-wise tree growth strategy.

---

← Model Building | Next: Interview Questions →