Linear regression

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.linear_model import LinearRegression

df = pd.read_csv(DATA / "uy_ijara.csv")
num = ["xonalar", "maydon_m2", "qavat", "qavatlar_soni", "metro_km", "mebel"]
X, y = df[num], df["narx_usd"]
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=42)

model = LinearRegression().fit(X_tr, y_tr)
print(f"bias (b) = {model.intercept_:.2f}")
print()
print("Coefficients:")
for name, k in zip(num, model.coef_):
    print(f"  {name:16s} {k:>9.3f}")

notes = []
for name, k in zip(num, model.coef_):
    direction = "goes up" if k > 0 else "goes down"
    unit = {"xonalar": "1 room", "maydon_m2": "1 m2", "qavat": "1 floor",
            "qavatlar_soni": "1 floor", "metro_km": "1 km", "mebel": "if furnished"}[name]
    notes.append({"Feature": name, "Coefficient": round(k, 2),
                  "How to read it": f"{unit} -> price {direction} by {abs(k):.1f} USD"})
show(pd.DataFrame(notes))

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

scaled = Pipeline([("sc", StandardScaler()), ("m", LinearRegression())]).fit(X_tr, y_tr)
k_scaled = scaled.named_steps["m"].coef_

compare = pd.DataFrame({
    "feature": num,
    "raw coefficient": model.coef_.round(2),
    "scaled": k_scaled.round(2),
}).assign(importance=lambda d: d.scaled.abs().round(2)).sort_values("importance", ascending=False)
show(compare)
print()
print("Now the correct conclusion: the most influential feature is", compare.iloc[0].feature)

from sklearn.metrics import r2_score, root_mean_squared_error, mean_absolute_error
from sklearn.dummy import DummyRegressor

pred = model.predict(X_te)
baseline = DummyRegressor(strategy="mean").fit(X_tr, y_tr)
b_pred = baseline.predict(X_te)

print(f"{'metric':<8} {'baseline (mean)':>17} {'model':>12} {'how to read it'}")
print("-" * 72)
print(f"{'R2':<8} {r2_score(y_te, b_pred):>17.4f} {r2_score(y_te, pred):>12.4f}   1.0 = perfect, 0 = same as the mean")
print(f"{'RMSE':<8} {root_mean_squared_error(y_te, b_pred):>17.1f} {root_mean_squared_error(y_te, pred):>12.1f}   USD, punishes big errors more")
print(f"{'MAE':<8} {mean_absolute_error(y_te, b_pred):>17.1f} {mean_absolute_error(y_te, pred):>12.1f}   USD, the average error")

show(pd.DataFrame([
    {"Metric": "R2", "What it says": "What share of the variance is explained",
     "Who to tell": "An engineer", "Weakness": "No business meaning"},
    {"Metric": "RMSE", "What it says": "The error, in the target's unit (USD)",
     "Who to tell": "An engineer", "Weakness": "Outliers pull it"},
    {"Metric": "MAE", "What it says": "The average error, in USD",
     "Who to tell": "THE MANAGER", "Weakness": "Goes easy on big errors"},
    {"Metric": "MAPE", "What it says": "The error as a percentage",
     "Who to tell": "The manager", "Weakness": "Explodes when the target is near zero"},
]))

mae = mean_absolute_error(y_te, pred)
print(f"The sentence for the manager:")
print(f'  "The model predicts rent with an average error of {mae:.0f} USD."')
print(f'  (the average rent is {y.mean():.0f} USD, so an error of ~{mae/y.mean():.1%})')

import matplotlib.pyplot as plt

resid = y_te - pred

figure, axes = plt.subplots(1, 3, figsize=(16, 4.4))

axes[0].scatter(pred, resid, alpha=.45, s=22, color="#1a73e8")
axes[0].axhline(0, color="red", ls="--", lw=1.5)
axes[0].set_xlabel("Prediction"); axes[0].set_ylabel("Residual (truth - prediction)")
axes[0].set_title("Residual vs prediction", fontweight="bold")

axes[1].hist(resid, bins=32, color="#43a047", alpha=.8)
axes[1].axvline(0, color="red", ls="--")
axes[1].set_xlabel("Residual"); axes[1].set_title("Distribution of residuals", fontweight="bold")

axes[2].scatter(y_te, pred, alpha=.45, s=22, color="#8e24aa")
limits = [min(y_te.min(), pred.min()), max(y_te.max(), pred.max())]
axes[2].plot(limits, limits, "r--", lw=1.5)
axes[2].set_xlabel("Truth"); axes[2].set_ylabel("Prediction")
axes[2].set_title("Truth vs prediction", fontweight="bold")

for a in axes: a.grid(alpha=.3)
fig("qoldiq.png", caption="Left: residuals should scatter randomly around zero. Middle: they should look roughly normal. Right: the points should be close to the diagonal.")

show(pd.DataFrame([
    {"What you see": "A random cloud around zero", "Diagnosis": "GOOD: the model has the right shape"},
    {"What you see": "A curve (a U or an arch)", "Diagnosis": "The relationship is NOT LINEAR -> a tree or a polynomial"},
    {"What you see": "A funnel (wider on the right)", "Diagnosis": "Heteroscedasticity -> try log(y)"},
    {"What you see": "A few points very far away", "Diagnosis": "Outliers -> read them, then decide"},
]))

# we check it with numbers
print(f"Mean of the residuals : {resid.mean():>9.4f}   (should be close to zero)")
print(f"Std of the residuals  : {resid.std():>9.2f}")
print()
# does the error grow as the prediction grows (the funnel test)
group = pd.qcut(pred, 4, labels=["cheapest", "cheap", "pricey", "priciest"])
check = pd.DataFrame({"group": group, "abs_error": np.abs(resid.to_numpy())})
print("Mean absolute error by prediction group:")
print(check.groupby("group", observed=True).abs_error.mean().round(1).to_string())

# add one strange value to the data
X_broken = X_tr.copy()
y_broken = y_tr.copy()
y_broken.iloc[0] = 15000.0          # a data-entry error: 15 000 USD/month

from sklearn.linear_model import HuberRegressor

rows = []
for name, m in [("LinearRegression (MSE)", LinearRegression()),
                ("HuberRegressor (close to MAE)", HuberRegressor(max_iter=2000))]:
    clean_m = Pipeline([("sc", StandardScaler()), ("m", m)]).fit(X_tr, y_tr)
    broken_m = Pipeline([("sc", StandardScaler()), ("m", m)]).fit(X_broken, y_broken)
    rows.append({
        "Model": name,
        "MAE (clean data)": round(mean_absolute_error(y_te, clean_m.predict(X_te)), 1),
        "MAE (with 1 strange value)": round(mean_absolute_error(y_te, broken_m.predict(X_te)), 1),
    })
show(pd.DataFrame(rows))

from sklearn.preprocessing import PolynomialFeatures

print(f"{'model':<34} {'CV R2':>10} {'±':>8}")
print("-" * 54)
for name, pipe in [
    ("Linear (6 features)", Pipeline([("sc", StandardScaler()), ("m", LinearRegression())])),
    ("Polynomial degree 2", Pipeline([("p", PolynomialFeatures(2, include_bias=False)),
                                      ("sc", StandardScaler()), ("m", LinearRegression())])),
    ("Polynomial degree 3", Pipeline([("p", PolynomialFeatures(3, include_bias=False)),
                                      ("sc", StandardScaler()), ("m", LinearRegression())])),
]:
    b = cross_val_score(pipe, X, y, cv=5, scoring="r2")
    print(f"{name:<34} {b.mean():>10.4f} {b.std():>8.4f}")

p2 = PolynomialFeatures(2, include_bias=False).fit(X)
p3 = PolynomialFeatures(3, include_bias=False).fit(X)
print(f"Number of features: 6 -> {p2.transform(X[:1]).shape[1]} (degree 2) -> {p3.transform(X[:1]).shape[1]} (degree 3)")
print()
print("The verdict, using the Chapter 10 rule: is the gap bigger than the fold spread?")

show(pd.DataFrame([
    {"Assumption": "The relationship is linear", "How you check it": "The residual plot: no curved pattern",
     "If it breaks": "A polynomial, a tree, or a log transform"},
    {"Assumption": "Residuals are independent", "How you check it": "A residual plot over time",
     "If it breaks": "Time-series methods"},
    {"Assumption": "Constant variance", "How you check it": "No funnel in the residual plot",
     "If it breaks": "log(y) or weighted regression"},
    {"Assumption": "Little multicollinearity", "How you check it": "The correlation matrix (Chapter 4)",
     "If it breaks": "Ridge (Chapter 17) or drop a column"},
    {"Assumption": "Few outliers", "How you check it": "boxplot",
     "If it breaks": "Huber or cleaning"},
]))

print("How unstable the coefficients are across different train splits:")
print(f"{'random_state':>13} {'xonalar':>10} {'maydon_m2':>12}")
print("-" * 38)
for rs in [0, 1, 7, 42, 99]:
    Xa, _, ya, _ = train_test_split(X, y, test_size=0.3, random_state=rs)
    m = Pipeline([("sc", StandardScaler()), ("m", LinearRegression())]).fit(Xa, ya)
    k = m.named_steps["m"].coef_
    print(f"{rs:>13} {k[0]:>10.2f} {k[1]:>12.2f}")

show(pd.DataFrame([
    {"Concept": "The equation", "In one sentence": "y = X @ w + b, the operation from Chapter 12"},
    {"Concept": "Reading a coefficient", "In one sentence": "'If it goes up by 1 unit, the target changes by this much'"},
    {"Concept": "Comparing coefficients", "In one sentence": "Scale first, or the conclusion comes out backwards"},
    {"Concept": "R2 / RMSE / MAE", "In one sentence": "Tell the manager the MAE: it's in USD"},
    {"Concept": "The residual plot", "In one sentence": "The most important diagnostic in regression"},
    {"Concept": "MSE and outliers", "In one sentence": "One broken row moves the model -> Huber"},
    {"Concept": "Polynomial features", "In one sentence": "A linear model can draw a curve"},
    {"Concept": "Multicollinearity", "In one sentence": "Unstable coefficients; the prediction stays stable"},
]))