Linear Regression — A Complete Worked Project

This lesson walks through a full regression project on a realistic dataset, applying every stage from Units 1–3: loading, EDA, cleaning, feature engineering, train-test split, modelling, evaluation, and diagnostics.

Problem statement: Predict a house's selling price from its characteristics, and identify which factors drive price most strongly.

Step 1 — Load and Inspect

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

np.random.seed(42)
n = 400

area      = np.random.normal(1500, 450, n).clip(500, 3500)
bedrooms  = np.clip(np.round(area / 550 + np.random.normal(0, 0.6, n)), 1, 6)
age       = np.random.randint(0, 40, n)
distance  = np.random.uniform(1, 25, n)          # km from city centre
location  = np.random.choice(["Prime", "Suburb", "Outskirts"], n, p=[0.3, 0.45, 0.25])

loc_premium = pd.Series(location).map({"Prime": 900000, "Suburb": 300000, "Outskirts": 0})
price = (2800 * area + 250000 * bedrooms - 22000 * age
         - 45000 * distance + loc_premium.values
         + np.random.normal(0, 250000, n))

df = pd.DataFrame({
    "area_sqft": area.round(0), "bedrooms": bedrooms, "age_years": age,
    "distance_km": distance.round(2), "location": location,
    "price": price.round(0).clip(500000, None),
})

print(df.shape)          # (400, 6)
print(df.head())
print(df.info())
print(df.describe().round(2))

Step 2 — Exploratory Data Analysis

print("Missing values:\n", df.isnull().sum())
print("Duplicates:", df.duplicated().sum())

numeric = df.select_dtypes(include=np.number)
print("\nCorrelation with price:")
print(numeric.corr()["price"].sort_values(ascending=False).round(3))
# price          1.000
# area_sqft      0.847
# bedrooms       0.611
# age_years     -0.278
# distance_km   -0.351

print("\nMean price by location:")
print(df.groupby("location")["price"].agg(["count", "mean", "median"]).round(0))

fig, axes = plt.subplots(2, 3, figsize=(17, 9))

axes[0, 0].hist(df["price"] / 1e5, bins=30, color="#168B99", edgecolor="white")
axes[0, 0].set_title(f"Price Distribution (skew = {df['price'].skew():.2f})")
axes[0, 0].set_xlabel("Price (Rs lakh)")

axes[0, 1].scatter(df["area_sqft"], df["price"] / 1e5, alpha=0.5, color="#10b981")
axes[0, 1].set_title("Price vs Area — strongest predictor")
axes[0, 1].set_xlabel("Area (sqft)"); axes[0, 1].set_ylabel("Price (Rs lakh)")

axes[0, 2].scatter(df["distance_km"], df["price"] / 1e5, alpha=0.5, color="#f59e0b")
axes[0, 2].set_title("Price vs Distance — negative relationship")

sns.boxplot(data=df, x="location", y=df["price"] / 1e5, ax=axes[1, 0], palette="Set2")
axes[1, 0].set_title("Price by Location")

sns.boxplot(data=df, x="bedrooms", y=df["price"] / 1e5, ax=axes[1, 1], palette="Set3")
axes[1, 1].set_title("Price by Bedrooms")

sns.heatmap(numeric.corr(), annot=True, cmap="coolwarm", center=0,
            fmt=".2f", ax=axes[1, 2], square=True)
axes[1, 2].set_title("Correlation Heatmap")

plt.tight_layout(); plt.show()

Step 3 — Feature Engineering and Encoding

model_df = df.copy()

# Derived features discovered during EDA
model_df["price_per_sqft"] = (model_df["price"] / model_df["area_sqft"]).round(2)
model_df["area_per_bedroom"] = (model_df["area_sqft"] / model_df["bedrooms"]).round(2)
model_df["is_new"] = (model_df["age_years"] <= 5).astype(int)

# One-hot encode the nominal 'location' column (drop_first avoids the dummy trap)
model_df = pd.get_dummies(model_df, columns=["location"], drop_first=True, dtype=int)

print(model_df.columns.tolist())
# ['area_sqft', 'bedrooms', 'age_years', 'distance_km', 'price',
#  'price_per_sqft', 'area_per_bedroom', 'is_new',
#  'location_Prime', 'location_Suburb']

# price_per_sqft is derived FROM the target — including it would leak the answer.
model_df = model_df.drop(columns=["price_per_sqft"])

Data leakage alert: price_per_sqft was computed from the target. Leaving it in would give a near-perfect R² in testing and a useless model in production, because at prediction time the price is exactly what you do not know.

Step 4 — Train-Test Split

X = model_df.drop(columns="price")
y = model_df["price"]

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)
print(f"Training: {X_train.shape}, Test: {X_test.shape}")
# Training: (320, 8), Test: (80, 8)

Step 5 — Train the Model

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

coeffs = pd.DataFrame({
    "feature": X.columns,
    "coefficient": model.coef_.round(2),
}).sort_values("coefficient", key=abs, ascending=False)

print("Intercept:", round(model.intercept_, 2))
print(coeffs.to_string(index=False))
# The 'location_Prime' coefficient is the price premium of a Prime location
# relative to the dropped baseline (Outskirts), holding everything else fixed.

Step 6 — Evaluate

y_train_pred = model.predict(X_train)
y_test_pred  = model.predict(X_test)

def report(y_true, y_pred, label):
    print(f"\n--- {label} ---")
    print(f"MAE  : Rs {mean_absolute_error(y_true, y_pred):,.0f}")
    print(f"RMSE : Rs {np.sqrt(mean_squared_error(y_true, y_pred)):,.0f}")
    print(f"R²   : {r2_score(y_true, y_pred):.4f}")
    mape = np.mean(np.abs((y_true - y_pred) / y_true)) * 100
    print(f"MAPE : {mape:.2f}%")

report(y_train, y_train_pred, "TRAINING SET")
report(y_test,  y_test_pred,  "TEST SET")

# Train and test R² being close together means the model GENERALISES —
# a large gap would signal overfitting.

# Cross-validation for a stable performance estimate
cv_scores = cross_val_score(LinearRegression(), X, y, cv=5, scoring="r2")
print("\n5-fold CV R²:", cv_scores.round(4))
print(f"Mean R² = {cv_scores.mean():.4f} (+/- {cv_scores.std():.4f})")

Step 7 — Diagnostics

residuals = y_test - y_test_pred

fig, axes = plt.subplots(1, 3, figsize=(17, 4.5))

axes[0].scatter(y_test / 1e5, y_test_pred / 1e5, alpha=0.6, color="#168B99")
lims = [min(y_test.min(), y_test_pred.min()) / 1e5,
        max(y_test.max(), y_test_pred.max()) / 1e5]
axes[0].plot(lims, lims, "r--", linewidth=2)
axes[0].set_xlabel("Actual price (Rs lakh)"); axes[0].set_ylabel("Predicted (Rs lakh)")
axes[0].set_title(f"Actual vs Predicted (R² = {r2_score(y_test, y_test_pred):.3f})")

axes[1].scatter(y_test_pred / 1e5, residuals / 1e5, alpha=0.6, color="#10b981")
axes[1].axhline(0, color="#ef4444", linestyle="--")
axes[1].set_xlabel("Fitted (Rs lakh)"); axes[1].set_ylabel("Residual (Rs lakh)")
axes[1].set_title("Residuals vs Fitted — want a random cloud")

axes[2].hist(residuals / 1e5, bins=25, color="#6366f1", edgecolor="white")
axes[2].set_title("Residual Distribution — want bell shaped")
axes[2].set_xlabel("Residual (Rs lakh)")

plt.tight_layout(); plt.show()

from scipy import stats
print("Shapiro-Wilk p on residuals:", round(stats.shapiro(residuals)[1], 4))

Step 8 — Predict for a New House

new_house = pd.DataFrame([{
    "area_sqft": 1800, "bedrooms": 3, "age_years": 4, "distance_km": 8.5,
    "area_per_bedroom": 1800 / 3, "is_new": 1,
    "location_Prime": 1, "location_Suburb": 0,
}])[X.columns]        # reorder columns to match training exactly

predicted = model.predict(new_house)[0]
print(f"Predicted price: Rs {predicted:,.0f}  (Rs {predicted/1e5:.2f} lakh)")

# A rough prediction interval from the residual spread
rmse = np.sqrt(mean_squared_error(y_test, y_test_pred))
print(f"Approx 95% interval: Rs {predicted - 1.96*rmse:,.0f} to Rs {predicted + 1.96*rmse:,.0f}")

Column order matters. scikit-learn matches features by position, not by name — reindexing new data to X.columns prevents silently wrong predictions.

Step 9 — Business Interpretation

insights = pd.DataFrame({
    "feature": X.columns,
    "coefficient": model.coef_,
}).assign(impact=lambda d: d["coefficient"].abs()).sort_values("impact", ascending=False)

print("\nKEY DRIVERS OF PRICE\n" + "=" * 55)
for _, row in insights.head(5).iterrows():
    direction = "increases" if row["coefficient"] > 0 else "decreases"
    print(f"{row['feature']:20} each +1 unit {direction} price by "
          f"Rs {abs(row['coefficient']):,.0f}")

Step 10 — Save and Reuse the Model

import joblib

joblib.dump(model, "house_price_model.pkl")
joblib.dump(list(X.columns), "model_features.pkl")

# Later, in production:
# loaded_model = joblib.load("house_price_model.pkl")
# feature_order = joblib.load("model_features.pkl")
# prediction = loaded_model.predict(new_data[feature_order])

Project Checklist

That completes supervised learning. The remaining Unit 3 lessons cover unsupervised techniques — clustering and association rule mining.