Practicals: Simple and Multiple Linear Regression

The theory is complete — now you run it on real data. The California Housing dataset has 20,640 census blocks, 8 features, and enough messiness (outliers, correlated features, heteroscedasticity) to expose everything a single-feature or toy dataset hides. The goal is concrete numbers you can verify and a residual analysis that tells you what the linear model still can't do.

Dataset Exploration

from sklearn.datasets import fetch_california_housing
import pandas as pd
import numpy as np

data = fetch_california_housing()
df = pd.DataFrame(data.data, columns=data.feature_names)
df['MedHouseVal'] = data.target

print(df.shape)
print(df.describe())
print(df.isnull().sum())

No null values. Two columns warrant attention:

• AveRooms: max = 141.9 (mean 5.4) — a handful of census blocks have extreme values, likely data errors.
• AveOccup: max = 1243.3 (mean 3.1) — severe outliers. Some sparsely populated blocks with a few homes show absurd occupancy ratios.

These outliers will inflate RMSE and create heteroscedasticity in residuals.

Simple Linear Regression — Single Feature

Start with MedInc (median income), which has the highest expected correlation with house value:

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score

X_simple = df[['MedInc']].values
y = df['MedHouseVal'].values

X_train, X_test, y_train, y_test = train_test_split(
    X_simple, y, test_size=0.2, random_state=42
)

model_simple = LinearRegression()
model_simple.fit(X_train, y_train)

print(f"Intercept:    {model_simple.intercept_:.4f}")
print(f"Coef[MedInc]: {model_simple.coef_[0]:.4f}")

For each unit increase in median income ($10k), house value increases by 0.4163 × $100k = $41,630.

y_pred_test = model_simple.predict(X_test)
print(f"Test RMSE: {np.sqrt(mean_squared_error(y_test, y_pred_test)):.4f}")
print(f"Test R²:   {r2_score(y_test, y_pred_test):.4f}")

$R^2=0.47$: income alone explains 47% of house value variance. RMSE of 0.83 in target units ($100k) means predictions are off by roughly $83k on average — a wide band.

The wide residual scatter shows that income captures the trend but leaves substantial unexplained variance — geography and house characteristics matter.

Multiple Linear Regression — All 8 Features

X_multi = df[data.feature_names].values

X_train, X_test, y_train, y_test = train_test_split(
    X_multi, y, test_size=0.2, random_state=42
)

model_multi = LinearRegression()
model_multi.fit(X_train, y_train)

print("Coefficients:")
for name, coef in zip(data.feature_names, model_multi.coef_):
    print(f"  {name:12s}: {coef:.4f}")
print(f"Intercept: {model_multi.intercept_:.4f}")

Notable signs:

• AveRooms: −0.107 but AveBedrms: +0.645 — these are highly correlated, and multicollinearity caused the sign flip. Don't interpret these individually.
• Latitude/Longitude: both negative — house values drop as you move north or east from coastal California. The model is encoding geography.
• Population: ≈0 — census block population, once other features are included, adds no marginal signal.

y_pred_multi = model_multi.predict(X_test)
print(f"Test RMSE: {np.sqrt(mean_squared_error(y_test, y_pred_multi)):.4f}")
print(f"Test R²:   {r2_score(y_test, y_pred_multi):.4f}")

$R^2$ improved from 0.47 (simple) to 0.60 (multiple) — 13 points from 7 additional features. RMSE dropped from 0.835 to 0.736 ($83.5k → $73.6k average error). The improvement is real but the ceiling for a linear model on this dataset is around $R^2=0.60$.

Residual Analysis

import matplotlib.pyplot as plt

residuals = y_test - y_pred_multi

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

plt.figure(figsize=(10, 4))

plt.subplot(1, 2, 1)
plt.hist(residuals, bins=50)
plt.xlabel('Residual')
plt.ylabel('Count')
plt.title('Residual Distribution')

plt.subplot(1, 2, 2)
plt.scatter(y_pred_multi, residuals, alpha=0.3, s=5)
plt.axhline(0, color='red')
plt.xlabel('Fitted Values')
plt.ylabel('Residuals')
plt.title('Residuals vs Fitted')
plt.tight_layout()
plt.show()

Two observations:

• Residual distribution: roughly symmetric around zero (mean ≈ 0 ✓) but with heavy left tail — high-value houses are underpredicted.
• Residuals vs Fitted: a fan shape for high fitted values ($>3.5$). This is heteroscedasticity — variance of errors grows with predicted value. One cause: the dataset caps house values at $500k ($5.0 in target units), so wealthy neighborhoods are systematically compressed.

Model Comparison

Related Concepts and Honest Limitations

The AveRooms/AveBedrms sign flip is multicollinearity in action. Their correlation is ~0.85 — the model assigns arbitrary individual signs to explain variance that belongs to both. Ridge regression (post 14) stabilizes this by shrinking both toward zero proportionally.

The 0.60 ceiling for $R^2$ on California Housing is a property of the dataset, not a flaw in linear regression. The true relationship between geography, income, and house price is non-linear — tree models (Random Forest, XGBoost) reach $R^2\approx 0.82$ on the same data. Linear regression is an honest baseline, not the end state.

Test Your Understanding

1. The multiple regression $R^2$ improved from 0.47 to 0.60 by adding 7 features. If you added 100 more randomly generated noise features, would $R^2$ increase, decrease, or stay the same? What about adjusted $R^2$?

2. AveRooms has a negative coefficient (−0.107) despite larger homes typically costing more. Explain how multicollinearity with AveBedrms produces this sign.

3. The mean residual is −0.0024 (near zero). Why is unbiasedness guaranteed by OLS on the training set? Is it guaranteed on the test set?

4. The residuals vs fitted plot shows a fan shape for high values. Name two transformations of $y$ that could reduce heteroscedasticity, and explain the trade-off of each.

5. Both AveRooms and AveBedrms have unreliable coefficients due to multicollinearity. If you had to drop one of these features, how would you decide which one to keep?