Linear Regression

Linear regression models the relationship between a dependent variable Y and one or more independent variables X by fitting a straight line (or hyperplane) that minimises the sum of squared errors.

Simple Linear Regression

   Ŷ = a + bX              or        Y = β₀ + β₁X + ε

   Ŷ = predicted value of Y
   a (β₀) = intercept — the value of Y when X = 0
   b (β₁) = slope — the change in Y for a one-unit increase in X
   ε = random error term

The Least Squares Method

The "best" line is the one minimising the sum of squared residuals:

   Minimise:  SSE = Σ (yᵢ − ŷᵢ)²  =  Σ (yᵢ − a − bxᵢ)²

Solving the normal equations gives:

           Σ(x − x̄)(y − ȳ)          n·Σxy − Σx·Σy            Cov(x, y)
   b  =  ───────────────────  =  ───────────────────────  =  ───────────
              Σ(x − x̄)²             n·Σx² − (Σx)²              Var(x)

   a  =  ȳ − b·x̄

Also useful:      b = r × (σy / σx)

The regression line always passes through the point (x̄, ȳ) — a handy check on your arithmetic.

Worked Example — Complete Hand Calculation

Hours studied (X) vs marks scored (Y):

Step 1 — Means

x̄ = 30/5 = 6          ȳ = 310/5 = 62

Step 2 — Slope

       n·ΣXY − ΣX·ΣY        5(2080) − (30)(310)      10400 − 9300      1100
  b = ─────────────────  =  ─────────────────────  =  ────────────  =  ──────  =  5.5
       n·ΣX² − (ΣX)²          5(220) − (30)²          1100 − 900        200

Step 3 — Intercept

  a = ȳ − b·x̄ = 62 − 5.5(6) = 62 − 33 = 29

Step 4 — Regression equation

  Ŷ = 29 + 5.5X

Interpretation:

• Slope 5.5 — each additional study hour is associated with 5.5 more marks
• Intercept 29 — a student who studies zero hours is predicted to score 29 marks

Step 5 — Prediction

For X = 7 hours:   Ŷ = 29 + 5.5(7) = 29 + 38.5 = 67.5 marks

Step 6 — Goodness of fit

       SS_res         20
R² = 1 − ────── = 1 − ────── = 1 − 0.01626 = 0.9837
       SS_tot        1230

=> 98.37% of the variation in marks is explained by study hours.

Check:  r = 0.9919 (computed in the correlation lesson),  r² = 0.9839  ✓
        (small difference is rounding)

RMSE = √(SS_res / n) = √(20/5) = √4 = 2.0 marks

Note: the residuals sum to zero — always true for a least-squares fit with an intercept.

Multiple Linear Regression

   Ŷ = β₀ + β₁X₁ + β₂X₂ + … + βₖXₖ

Each βᵢ is the change in Y per one-unit change in Xᵢ, HOLDING ALL OTHER
PREDICTORS CONSTANT — the "partial" or "ceteris paribus" interpretation.

Matrix solution (Ordinary Least Squares):

   β = (XᵀX)⁻¹ Xᵀy

Python

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_squared_error

X = np.array([2, 4, 6, 8, 10]).reshape(-1, 1)
y = np.array([40, 50, 65, 70, 85])

model = LinearRegression()
model.fit(X, y)

print(f"Slope (b)     = {model.coef_[0]:.4f}")        # 5.5000
print(f"Intercept (a) = {model.intercept_:.4f}")      # 29.0000
print(f"Equation: Y = {model.intercept_:.2f} + {model.coef_[0]:.2f}X")

y_pred = model.predict(X)
print(f"R²   = {r2_score(y, y_pred):.4f}")                        # 0.9837
print(f"RMSE = {np.sqrt(mean_squared_error(y, y_pred)):.4f}")     # 2.0000

print(f"Prediction for 7 hours: {model.predict([[7]])[0]:.2f}")   # 67.50

# Computing the coefficients manually — exactly the exam formula
x = np.array([2, 4, 6, 8, 10], dtype=float)
y = np.array([40, 50, 65, 70, 85], dtype=float)
n = len(x)

b = (n * (x * y).sum() - x.sum() * y.sum()) / (n * (x ** 2).sum() - x.sum() ** 2)
a = y.mean() - b * x.mean()
print(f"Manual: b = {b:.4f}, a = {a:.4f}")     # b = 5.5000, a = 29.0000

# Alternative: b = r × (σy / σx)
r = np.corrcoef(x, y)[0, 1]
b_alt = r * (y.std(ddof=1) / x.std(ddof=1))
print(f"Via correlation: b = {b_alt:.4f}")      # 5.5000  ✓

# Visualising the fit and the residuals
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

axes[0].scatter(X, y, s=90, color="#168B99", zorder=3, label="Actual")
axes[0].plot(X, y_pred, color="#ef4444", linewidth=2,
             label=f"Ŷ = {model.intercept_:.1f} + {model.coef_[0]:.1f}X")
for xi, yi, ypi in zip(X.ravel(), y, y_pred):        # draw the residuals
    axes[0].plot([xi, xi], [yi, ypi], color="gray", linestyle=":", zorder=2)
axes[0].scatter([x.mean()], [y.mean()], marker="X", s=200, color="#f59e0b",
                zorder=4, label="(x̄, ȳ) — line passes through this")
axes[0].set_xlabel("Hours Studied"); axes[0].set_ylabel("Marks")
axes[0].set_title(f"Least Squares Fit (R² = {r2_score(y, y_pred):.4f})")
axes[0].legend(); axes[0].grid(alpha=0.3)

axes[1].scatter(y_pred, y - y_pred, s=90, color="#10b981")
axes[1].axhline(0, color="#ef4444", linestyle="--")
axes[1].set_xlabel("Fitted values"); axes[1].set_ylabel("Residuals")
axes[1].set_title("Residual Plot")
axes[1].grid(alpha=0.3)

plt.tight_layout(); plt.show()

# MULTIPLE linear regression
df = pd.DataFrame({
    "hours":      [2, 4, 6, 8, 10, 3, 7, 9, 5, 11],
    "attendance": [60, 70, 80, 85, 95, 65, 82, 90, 75, 97],
    "prev_score": [55, 60, 68, 72, 80, 58, 70, 78, 65, 85],
    "marks":      [40, 50, 65, 70, 85, 45, 68, 80, 58, 90],
})

X = df[["hours", "attendance", "prev_score"]]
y = df["marks"]

mlr = LinearRegression().fit(X, y)

print("Intercept:", round(mlr.intercept_, 4))
print(pd.DataFrame({"feature": X.columns, "coefficient": mlr.coef_.round(4)}))
print(f"R² = {mlr.score(X, y):.4f}")

# Adjusted R² — the honest metric when comparing models with different k
n, k = X.shape
r2 = mlr.score(X, y)
adj_r2 = 1 - (1 - r2) * (n - 1) / (n - k - 1)
print(f"Adjusted R² = {adj_r2:.4f}")

# statsmodels gives the full statistical summary that exams reference
import statsmodels.api as sm

X_sm = sm.add_constant(df[["hours", "attendance", "prev_score"]])
ols = sm.OLS(df["marks"], X_sm).fit()
print(ols.summary())
# Reports: coefficients, standard errors, t-statistics, p-values for each
# predictor, R², adjusted R², F-statistic, confidence intervals, and
# diagnostics (Durbin-Watson, Jarque-Bera).
# A predictor with p > 0.05 is not making a statistically significant
# contribution once the others are accounted for.

Interpreting the Coefficients

Extrapolation Warning

The model is only valid over the range of X in the training data. Our equation predicts Ŷ = 29 + 5.5(30) = 194 marks for 30 study hours — impossible on a 100-mark exam. Never extrapolate far beyond the observed range.

Regularised Regression — When Predictors Are Many or Correlated

from sklearn.linear_model import Ridge, Lasso
from sklearn.preprocessing import StandardScaler

X_scaled = StandardScaler().fit_transform(X)

ridge = Ridge(alpha=1.0).fit(X_scaled, y)
lasso = Lasso(alpha=1.0).fit(X_scaled, y)

print(pd.DataFrame({
    "feature": X.columns,
    "OLS":   LinearRegression().fit(X_scaled, y).coef_.round(3),
    "Ridge": ridge.coef_.round(3),
    "Lasso": lasso.coef_.round(3),
}))
# Lasso drives weak predictors' coefficients to exactly 0.0

The next lesson applies everything here to a realistic end-to-end regression project.