Previous lesson. Advanced Regression Topics · Next lesson. Regression Quick Reference Guide

Complete Regression Project Tutorial

Walk a house-price regression from load to tuned Ridge. You will treat scaling, leakage, and error metrics the way a production ML review would.

Project Overview
Why this matters in production
Concept to application
Stand-out signal
Step 1: Data Loading and Exploration
Step 2: Data Cleaning and Preprocessing
Step 3: Feature Engineering
Step 4: Model Training
Step 5: Model Evaluation
Step 6: Model Improvement
Step 7: Final Model and Predictions

Project Overview

Project: Predict House Prices

Dataset: California Housing Dataset (built into scikit-learn)

Goal: Build a regression model to predict median house values

Type: Multiple Linear Regression with Regularization

Difficulty: Intermediate

Time: 1-2 hours

Why this matters in production

House price models look friendly until someone ships them. Then you meet target capping, spatial leakage, and stakeholders who hear “R² = 0.6” as “60% accurate.” Before you tune alpha, decide the decision the number supports (list price band, risk flag, ranking) and what a costly miss looks like in dollars.

Concept to application

Ridge shrinks coefficients so correlated features do not explode. That is the same idea behind many tabular pricing and demand models: linear structure you can audit, plus regularization so the fit survives noisy columns. Scaling before Ridge is not optional trivia. Without it, alpha does not mean what you think across feature ranges.

Stand-out signal

Show a short model card in your write-up: train/test protocol, whether geography leaked, RMSE in dollars, and one failure case (for example, luxury tails). Interviewers remember that more than another GridSearch screenshot. Outcomes still vary by role and market. This project builds evidence. It does not guarantee a job.

Step 1: Data Loading and Exploration

Load Data

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
from sklearn.model_selection import cross_val_score, GridSearchCV
import warnings
warnings.filterwarnings('ignore')

# Load California Housing dataset
housing = fetch_california_housing()
df = pd.DataFrame(housing.data, columns=housing.feature_names)
df['MedHouseVal'] = housing.target

print("Dataset loaded successfully!")
print(f"Shape: {df.shape}")
print(f"\nFirst few rows:")
print(df.head())

Basic Statistics

print("Dataset Info:")
print(df.info())

print("\nSummary Statistics:")
print(df.describe())

print("\nMissing Values:")
print(df.isnull().sum())

print("\nTarget Variable (MedHouseVal) Statistics:")
print(df['MedHouseVal'].describe())

Visualizations

# Distribution of target variable
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.hist(df['MedHouseVal'], bins=50, edgecolor='black')
plt.xlabel('Median House Value')
plt.ylabel('Frequency')
plt.title('Distribution of House Values')

plt.subplot(1, 2, 2)
plt.boxplot(df['MedHouseVal'])
plt.ylabel('Median House Value')
plt.title('Box Plot of House Values')

plt.tight_layout()
plt.show()

# Correlation heatmap
plt.figure(figsize=(10, 8))
correlation = df.corr()
sns.heatmap(correlation, annot=True, cmap='coolwarm', center=0, 
            square=True, fmt='.2f')
plt.title('Feature Correlation Matrix')
plt.tight_layout()
plt.show()

# Scatter plots of features vs target
fig, axes = plt.subplots(2, 4, figsize=(16, 8))
axes = axes.flatten()

for idx, feature in enumerate(housing.feature_names):
    axes[idx].scatter(df[feature], df['MedHouseVal'], alpha=0.3)
    axes[idx].set_xlabel(feature)
    axes[idx].set_ylabel('MedHouseVal')
    axes[idx].set_title(f'{feature} vs House Value')

plt.tight_layout()
plt.show()

Insights:

Check for outliers in target variable
Identify highly correlated features
Understand feature distributions
Check for non-linear relationships

Step 2: Data Cleaning and Preprocessing

Handle Missing Values

# Check for missing values
print("Missing values per column:")
print(df.isnull().sum())

# If there are missing values, handle them
# Option 1: Drop rows with missing values
# df = df.dropna()

# Option 2: Fill with median (for numerical)
# df = df.fillna(df.median())

# Option 3: Fill with mean
# df = df.fillna(df.mean())

Handle Outliers

Blindly dropping tails can erase the rare homes you care about most. Use IQR as a diagnostic first. Only filter when the business problem says those rows are sensor noise, not signal.

import pandas as pd


def remove_outliers_iqr(df, column):
    """Return rows whose column value lies inside the 1.5 IQR fences."""
    Q1 = df[column].quantile(0.25)
    Q3 = df[column].quantile(0.75)
    IQR = Q3 - Q1
    lower_bound = Q1 - 1.5 * IQR
    upper_bound = Q3 + 1.5 * IQR
    return df[(df[column] >= lower_bound) & (df[column] <= upper_bound)]

# Detect outliers using IQR method (optional cleaning step)
print(f"Original shape: {df.shape}")
df_clean = remove_outliers_iqr(df, "MedHouseVal")
print(f"After removing outliers: {df_clean.shape}")
print(
    f"Removed {len(df) - len(df_clean)} outliers "
    f"({100 * (len(df) - len(df_clean)) / len(df):.1f}%)"
)

# Visualize before and after
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].boxplot(df["MedHouseVal"])
axes[0].set_title("Before Outlier Removal")
axes[0].set_ylabel("MedHouseVal")

axes[1].boxplot(df_clean["MedHouseVal"])
axes[1].set_title("After Outlier Removal")
axes[1].set_ylabel("MedHouseVal")

plt.tight_layout()
plt.show()

df = df_clean  # Use cleaned data only if you accept the bias this introduces

Prepare Features and Target

# Separate features and target
X = df[housing.feature_names]
y = df['MedHouseVal']

print(f"Features shape: {X.shape}")
print(f"Target shape: {y.shape}")

Step 3: Feature Engineering

Check for Multicollinearity

# Calculate VIF (Variance Inflation Factor)
from statsmodels.stats.outliers_influence import variance_inflation_factor

def calculate_vif(X):
    vif_data = pd.DataFrame()
    vif_data["Feature"] = X.columns
    vif_data["VIF"] = [variance_inflation_factor(X.values, i) 
                       for i in range(X.shape[1])]
    return vif_data.sort_values('VIF', ascending=False)

vif_df = calculate_vif(X)
print("Variance Inflation Factors:")
print(vif_df)

# Features with VIF > 10 have multicollinearity
high_vif = vif_df[vif_df["VIF"] > 10]
if len(high_vif) > 0:
    print(f"\nWarning: {len(high_vif)} features have high VIF (>10)")
    print("Consider using regularization (Ridge/Lasso)")

Create New Features (Optional)

# Example: Create interaction features
# X['MedInc_x_AveRooms'] = X['MedInc'] * X['AveRooms']
# X['Population_x_HouseAge'] = X['Population'] * X['HouseAge']

# Example: Create polynomial features for specific features
# from sklearn.preprocessing import PolynomialFeatures
# poly = PolynomialFeatures(degree=2, include_bias=False, interacti>
# X_poly = poly.fit_transform(X[['MedInc', 'AveRooms']])

Split Data

# Split into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

print(f"Training set: {X_train.shape[0]} samples")
print(f"Test set: {X_test.shape[0]} samples")

Scale Features

# Scale features (important for regularization)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# Convert back to DataFrame for easier handling
X_train_scaled = pd.DataFrame(X_train_scaled, columns=housing.feature_names)
X_test_scaled = pd.DataFrame(X_test_scaled, columns=housing.feature_names)

print("Features scaled successfully!")

Step 4: Model Training

Baseline Model (Linear Regression)

# Train baseline model
baseline_model = LinearRegression()
baseline_model.fit(X_train_scaled, y_train)

# Predictions
y_train_pred = baseline_model.predict(X_train_scaled)
y_test_pred = baseline_model.predict(X_test_scaled)

# Evaluate
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("Baseline Model (Linear Regression):")
print(f"  Training RMSE: {train_rmse:.3f}")
print(f"  Test RMSE: {test_rmse:.3f}")
print(f"  Training R²: {train_r2:.3f}")
print(f"  Test R²: {test_r2:.3f}")

# Check for overfitting
if train_rmse < test_rmse * 0.9:
    print("  Warning: Possible overfitting (large gap between train and test)")

Ridge Regression

# Train Ridge regression
ridge_model = Ridge(alpha=1.0)
ridge_model.fit(X_train_scaled, y_train)

y_test_pred_ridge = ridge_model.predict(X_test_scaled)
test_rmse_ridge = np.sqrt(mean_squared_error(y_test, y_test_pred_ridge))
test_r2_ridge = r2_score(y_test, y_test_pred_ridge)

print("\nRidge Regression:")
print(f"  Test RMSE: {test_rmse_ridge:.3f}")
print(f"  Test R²: {test_r2_ridge:.3f}")

Lasso Regression

# Train Lasso regression
lasso_model = Lasso(alpha=0.1, max_iter=10000)
lasso_model.fit(X_train_scaled, y_train)

y_test_pred_lasso = lasso_model.predict(X_test_scaled)
test_rmse_lasso = np.sqrt(mean_squared_error(y_test, y_test_pred_lasso))
test_r2_lasso = r2_score(y_test, y_test_pred_lasso)

print("\nLasso Regression:")
print(f"  Test RMSE: {test_rmse_lasso:.3f}")
print(f"  Test R²: {test_r2_lasso:.3f}")

# Check feature selection
n > 0.001)
print(f"  Features used: {non_zero_coef} out of {len(housing.feature_names)}")

Step 5: Model Evaluation

Compare All Models

models = {
    'Linear Regression': (baseline_model, y_test_pred),
    'Ridge': (ridge_model, y_test_pred_ridge),
    'Lasso': (lasso_model, y_test_pred_lasso)
}

results = []
for name, (model, y_pred) in models.items():
    rmse = np.sqrt(mean_squared_error(y_test, y_pred))
    mae = mean_absolute_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    results.append({
        'Model': name,
        'RMSE': rmse,
        'MAE': mae,
        'R²': r2
    })

results_df = pd.DataFrame(results)
print("\nModel Comparison:")
print(results_df.to_string(index=False))

Residual Analysis

# Choose best model (lowest RMSE)
best_model_name = results_df.loc[results_df['RMSE'].idxmin(), 'Model']
best_model = models[best_model_name][0]
best_predicti>1]

print(f"\nBest Model: {best_model_name}")

# Residual analysis
residuals = y_test - best_predictions

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. Residuals vs Predicted
axes[0, 0].scatter(best_predictions, residuals, alpha=0.5)
axes[0, 0].axhline(y=0, color='r', linestyle='--')
axes[0, 0].set_xlabel('Predicted Values')
axes[0, 0].set_ylabel('Residuals')
axes[0, 0].set_title('Residuals vs Predicted')
axes[0, 0].grid(True)

# 2. Q-Q Plot
from scipy import stats
stats.probplot(residuals, dist="norm", plot=axes[0, 1])
axes[0, 1].set_title('Q-Q Plot (Normality Check)')
axes[0, 1].grid(True)

# 3. Histogram of Residuals
axes[1, 0].hist(residuals, bins=30, edgecolor='black')
axes[1, 0].set_xlabel('Residuals')
axes[1, 0].set_ylabel('Frequency')
axes[1, 0].set_title('Distribution of Residuals')
axes[1, 0].grid(True)

# 4. Actual vs Predicted
axes[1, 1].scatter(y_test, best_predictions, alpha=0.5)
axes[1, 1].plot([y_test.min(), y_test.max()], 
                [y_test.min(), y_test.max()], 'r--', linewidth=2)
axes[1, 1].set_xlabel('Actual Values')
axes[1, 1].set_ylabel('Predicted Values')
axes[1, 1].set_title('Actual vs Predicted')
axes[1, 1].grid(True)

plt.tight_layout()
plt.show()

# Residual statistics
print("\nResidual Statistics:")
print(f"  Mean: {residuals.mean():.6f} (should be ~0)")
print(f"  Std: {residuals.std():.3f}")
print(f"  Min: {residuals.min():.3f}")
print(f"  Max: {residuals.max():.3f}")

Feature Importance

# Feature importance (coefficients)
feature_importance = pd.DataFrame({
    'Feature': housing.feature_names,
    'Coefficient': best_model.coef_,
    'Abs_Coefficient': np.abs(best_model.coef_)
}).sort_values('Abs_Coefficient', ascending=False)

print("\nFeature Importance (sorted by absolute coefficient):")
print(feature_importance)

# Visualize
plt.figure(figsize=(10, 6))
plt.barh(feature_importance['Feature'], feature_importance['Coefficient'])
plt.xlabel('Coefficient Value')
plt.title('Feature Coefficients')
plt.tight_layout()
plt.show()

Step 6: Model Improvement

Hyperparameter Tuning

# Tune Ridge regression
ridge_params = {
    'alpha': [0.001, 0.01, 0.1, 1.0, 10.0, 100.0]
}

ridge_grid = GridSearchCV(
    Ridge(),
    ridge_params,
    cv=5,
    scoring='neg_mean_squared_error',
    return_train_score=True
)

ridge_grid.fit(X_train_scaled, y_train)

print("Ridge Regression - Best Parameters:")
print(f"  Alpha: {ridge_grid.best_params_['alpha']}")
print(f"  Best CV RMSE: {np.sqrt(-ridge_grid.best_score_):.3f}")

# Tune Lasso regression
lasso_params = {
    'alpha': [0.001, 0.01, 0.1, 1.0, 10.0]
}

lasso_grid = GridSearchCV(
    Lasso(max_iter=10000),
    lasso_params,
    cv=5,
    scoring='neg_mean_squared_error'
)

lasso_grid.fit(X_train_scaled, y_train)

print("\nLasso Regression - Best Parameters:")
print(f"  Alpha: {lasso_grid.best_params_['alpha']}")
print(f"  Best CV RMSE: {np.sqrt(-lasso_grid.best_score_):.3f}")

# Use best model
best_tuned_model = ridge_grid.best_estimator_
y_test_pred_tuned = best_tuned_model.predict(X_test_scaled)

test_rmse_tuned = np.sqrt(mean_squared_error(y_test, y_test_pred_tuned))
test_r2_tuned = r2_score(y_test, y_test_pred_tuned)

print(f"\nTuned Model Performance:")
print(f"  Test RMSE: {test_rmse_tuned:.3f}")
print(f"  Test R²: {test_r2_tuned:.3f}")

Cross-Validation

# Cross-validation scores
cv_scores = cross_val_score(
    best_tuned_model, X_train_scaled, y_train,
    cv=5, scoring='neg_mean_squared_error'
)

print(f"\nCross-Validation Results:")
print(f"  Mean RMSE: {np.sqrt(-cv_scores.mean()):.3f}")
rmse_folds = np.sqrt(-cv_scores)
print(f"  Std RMSE: {rmse_folds.std():.3f}")
print(
    f"  95% Confidence Interval: "
    f"[{rmse_folds.mean() - 1.96 * rmse_folds.std():.3f}, "
    f"{rmse_folds.mean() + 1.96 * rmse_folds.std():.3f}]"
)

Step 7: Final Model and Predictions

Final Evaluation

# Final model evaluation
def evaluate_model(y_true, y_pred, model_name):
    """Comprehensive model evaluation"""
    rmse = np.sqrt(mean_squared_error(y_true, y_pred))
    mae = mean_absolute_error(y_true, y_pred)
    r2 = r2_score(y_true, y_pred)
    
    print(f"\n{model_name} - Final Evaluation:")
    print(f"  RMSE: {rmse:.3f}")
    print(f"  MAE:  {mae:.3f}")
    print(f"  R²:   {r2:.3f}")
    
    # Interpretation
    print(f"\n  Interpretation:")
    print(f"    - Average prediction error: ${rmse*100000:.2f}")
    print(f"    - Model explains {r2*100:.1f}% of variance")
    
    return {'RMSE': rmse, 'MAE': mae, 'R²': r2}

final_metrics = evaluate_model(y_test, y_test_pred_tuned, "Final Model")

Make Predictions on New Data

# Example: Predict for new house
new_house = {
    'MedInc': 8.0,           # Median income
    'HouseAge': 20.0,        # House age
    'AveRooms': 5.0,         # Average rooms
    'AveBedrms': 1.0,        # Average bedrooms
    'Population': 2000.0,    # Population
    'AveOccup': 3.0,         # Average occupancy
    'Latitude': 34.0,        # Latitude
    'Longitude': -118.0      # Longitude
}

# Convert to DataFrame
new_house_df = pd.DataFrame([new_house])

# Scale features
new_house_scaled = scaler.transform(new_house_df)

# Predict
predicted_value = best_tuned_model.predict(new_house_scaled)[0]

print(f"\nPrediction for New House:")
for key, value in new_house.items():
    print(f"  {key}: {value}")
print(f"\nPredicted Median House Value: ${predicted_value*100000:,.2f}")

Save Model

import joblib

# Save model and scaler
joblib.dump(best_tuned_model, 'house_price_model.pkl')
joblib.dump(scaler, 'scaler.pkl')

print("\nModel and scaler saved successfully!")

# Load model (for future use)
# loaded_model = joblib.load('house_price_model.pkl')
# loaded_scaler = joblib.load('scaler.pkl')

Complete Code Summary

Self-contained pipeline used by the private hub snippet suite. Fit the scaler on train only. Tune Ridge with a small alpha grid.

# Complete Regression Project Pipeline
import numpy as np
import pandas as pd
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error, r2_score
import joblib

# 1. Load data
housing = fetch_california_housing()
df = pd.DataFrame(housing.data, columns=housing.feature_names)
df["MedHouseVal"] = housing.target

# 2. Split before fitting any transformer
X = df[housing.feature_names]
y = df["MedHouseVal"]
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# 3. Scale features (fit on train only)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# 4. Train and tune Ridge
ridge_params = {"alpha": [0.001, 0.01, 0.1, 1.0, 10.0, 100.0]}
ridge_grid = GridSearchCV(
    Ridge(), ridge_params, cv=5, scoring="neg_mean_squared_error"
)
ridge_grid.fit(X_train_scaled, y_train)

# 5. Evaluate on held-out test
best_model = ridge_grid.best_estimator_
y_pred = best_model.predict(X_test_scaled)
rmse = float(np.sqrt(mean_squared_error(y_test, y_pred)))
r2 = float(r2_score(y_test, y_pred))

print(f"RMSE: {rmse:.3f}")
print(f"R2: {r2:.3f}")
print(f"Best alpha: {ridge_grid.best_params_['alpha']}")

# 6. Persist artifacts (safe to run in a temp folder)
joblib.dump(best_model, "model.pkl")
joblib.dump(scaler, "scaler.pkl")

Key Takeaways

Explore before you fit - distributions and leakage beat fancy models
Treat outlier filters as a product choice - not a default delete
Scale before Ridge/Lasso - otherwise regularization is distorted
Report RMSE in dollars - stakeholders do not think in raw target units
Tune with cross-validation - one test score is not a design review
Save the scaler with the model - inference must match training transforms

Try next: Rebuild this with a residual plot and a short model card, then continue to 04-supervised-learning-classification.

All lessons in this module