{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# California housing: where the score comes from\n",
"\n",
"How well can a model predict the median house value of a California district from census data — and how much of the score belongs to the model, and how much to the way it is tested?\n",
"\n",
"*From [Pattern is Everything](https://patterniseverything.com/cases/housing-regression/). Data: California Housing — scikit-learn’s fetch_california_housing; originally StatLib, via Luís Torgo.*"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"# Fetch the data file if it is not next to this notebook (in JupyterLite it already is)\n",
"import os, urllib.request\n",
"if not os.path.exists('housing.csv'):\n",
" urllib.request.urlretrieve('https://patterniseverything.com/cases/data/housing.csv', 'housing.csv')"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Look before modelling\n",
"\n",
"Load the data and check the ranges. Three columns stop dead at a limit, the sign of values cut off (top-coded) when the data were published."
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"df = pd.read_csv('housing.csv')\n",
"df['value'] = df.median_house_value / 100_000 # in $100,000s\n",
"rooms = df.total_rooms / df.households\n",
"\n",
"n_rows = len(df)\n",
"at_value_cap = (df.median_house_value >= 500_001).sum() # top-coded at $500,001\n",
"at_age_cap = (df.housing_median_age == 52).sum() # top-coded at 52 years\n",
"at_income_cap = (df.median_income >= 15).sum() # top-coded at 15 ($150,000)\n",
"print(n_rows, 'districts')\n",
"print('at the value cap:', at_value_cap, f'({at_value_cap / n_rows:.1%})')\n",
"print('at the age cap:', at_age_cap, '| at the income cap:', at_income_cap)\n",
"print('rooms per household: median', round(rooms.median(), 1), '| max', round(rooms.max(), 1))"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"965 districts (4.7%) sit exactly at $500,001. Their true median value is somewhere above it, and no model trained on this column can learn by how much. The age and income columns are capped too, and a few districts average over 100 rooms per household — resort towns, where most houses are holiday homes with few permanent households."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. A baseline, then a test that matches the question\n",
"\n",
"Predicting the average value for every district misses by about $91k on average. Two quick models: a linear regression on income and housing features, and k-nearest neighbours that sees nothing but the location — it predicts the average value of the ten nearest districts. Each is scored twice: with random 5-fold cross-validation, and with folds made of whole regions (1°×1° squares of latitude and longitude, 55 of them), so the test districts lie in areas the model has not seen."
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"from sklearn.model_selection import KFold, PredefinedSplit, cross_val_score\n",
"from sklearn.linear_model import LinearRegression\n",
"from sklearn.neighbors import KNeighborsRegressor\n",
"\n",
"X = pd.DataFrame({\n",
" 'income': df.median_income, 'age': df.housing_median_age, 'rooms': rooms,\n",
" 'bedroom_share': df.total_bedrooms / df.total_rooms,\n",
" 'people_per_household': df.population / df.households,\n",
" 'lat': df.latitude, 'lon': df.longitude})\n",
"jitter = np.random.default_rng(0).normal(0, 1e-6, (len(X), 2))\n",
"X[['lat', 'lon']] += jitter # districts that share coordinates: makes \"nearest\" unique\n",
"y = df.value\n",
"baseline_mae = (y - y.mean()).abs().mean()\n",
"\n",
"block = np.floor(df.latitude).astype(int) * 1000 + np.floor(df.longitude).astype(int)\n",
"sizes = block.value_counts().sort_index().sort_values(ascending=False, kind='stable')\n",
"n_blocks = len(sizes)\n",
"load, fold = np.zeros(5), {}\n",
"for b, n in sizes.items(): # largest region first, each to the emptiest fold\n",
" k = int(load.argmin()); fold[b] = k; load[k] += n\n",
"random_cv = KFold(5, shuffle=True, random_state=0)\n",
"region_cv = PredefinedSplit(block.map(fold))\n",
"\n",
"def r2(model, cols):\n",
" a = cross_val_score(model, X[cols], y, cv=random_cv, scoring='r2').mean()\n",
" b = cross_val_score(model, X[cols], y, cv=region_cv, scoring='r2').mean()\n",
" return round(a, 3), round(b, 3)\n",
"\n",
"features = ['income', 'age', 'rooms', 'bedroom_share', 'people_per_household']\n",
"linear = r2(LinearRegression(), features)\n",
"knn = r2(KNeighborsRegressor(10), ['lat', 'lon'])\n",
"print('baseline MAE:', round(baseline_mae, 3))\n",
"print('R² (random folds, new regions)')\n",
"print(' linear, no location:', linear)\n",
"print(' 10 nearest districts:', knn)"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The linear model scores about the same either way (R² 0.55 and 0.49). The location-only model looks better than it — R² 0.79 with random folds — and collapses to 0.17 on new regions. With random folds, the nearest neighbours of a test district are its own neighbours in the training data; the model is not predicting, it is looking up the street. Nearby districts have similar prices (spatial autocorrelation), so random folds leak information whenever location is a feature."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. A stronger model, the same gap\n",
"\n",
"Gradient-boosted trees on all the features, location included — the kind of model that wins tabular benchmarks. It is fitted ten times, once per fold of each test."
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"from sklearn.ensemble import HistGradientBoostingRegressor\n",
"\n",
"boost = r2(HistGradientBoostingRegressor(random_state=0), list(X.columns))\n",
"print('gradient boosting, all features:', boost)"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"R² 0.84 with random folds, 0.68 on new regions. Both numbers are correct; they answer different questions. If the model will value houses in districts like the ones it was trained on — filling gaps on a map it already covers — the random score is the right one. If it will be used somewhere new, only the regional score says what to expect, and it is 19% lower."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4. What the cap does to the errors\n",
"\n",
"Hold out a random fifth of the districts, fit the boosted model on the rest, and look at the districts whose value is capped."
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"from sklearn.model_selection import train_test_split\n",
"\n",
"X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=0)\n",
"model = HistGradientBoostingRegressor(random_state=0).fit(X_tr, y_tr)\n",
"pred = pd.Series(model.predict(X_te), index=y_te.index)\n",
"err = (pred - y_te).abs()\n",
"\n",
"capped = y_te >= 5.00001\n",
"worst = err >= err.quantile(0.95) # the 5% largest errors\n",
"print('test MAE:', round(err.mean(), 3), '| capped districts in test:', capped.sum())\n",
"print('mean prediction for capped districts:', round(pred[capped].mean(), 2))\n",
"print('share of capped districts among the worst 5%:',\n",
" round((capped & worst).sum() / worst.sum(), 3))\n",
"print('highest prediction:', round(pred.max(), 2))"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The capped districts are 4% of the test set but 17% of the worst 5% of errors. For them the model predicts $445k on average, against a recorded $500,001 that is itself an underestimate. A metric averaged over all districts hides this: the error is concentrated at the expensive end, where a valuation mistake costs the most."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python"
}
},
"nbformat": 4,
"nbformat_minor": 5
}