Linear Regression with Multiple Regression (Expert) — Interactive Quiz

Definitions

Every key term in one or two sentences. Type to filter.

Term	Definition
Multiple linear regression	A model that predicts a numeric target as an intercept plus a weighted sum of two or more features, y = Xβ + ε, with weights estimated from data.
Design matrix (X)	The n × (p+1) matrix of feature values, one row per observation, usually with a leading column of ones for the intercept.
Intercept column	The column of ones in X; its coefficient β₀ shifts all predictions by a constant.
Coefficient (β_j)	The expected change in y per one-unit increase in feature x_j, holding all other features in the model constant.
Error term (ε)	The unobserved difference between the true y and the population regression function Xβ.
Residual (e)	The observed difference y − ŷ between actual and fitted values; the sample counterpart of ε.
Ordinary least squares (OLS)	The estimator that chooses β to minimise ‖y − Xβ‖², the sum of squared residuals.
Normal equations	The linear system XᵀXβ̂ = Xᵀy obtained by setting the gradient of the squared-error loss to zero.
Full column rank	The property that no column of X is a linear combination of the others; required for a unique OLS solution.
Hat matrix (H)	H = X(XᵀX)⁻¹Xᵀ, the projection matrix that maps y to fitted values ŷ = Hy.
Projection	The closest point in a subspace to a given vector; OLS fitted values are the orthogonal projection of y onto the column space of X.
Column space	The set of all vectors that can be written as Xβ for some β.
Leverage (hᵢᵢ)	The i-th diagonal entry of H; measures how unusual row i's feature values are and how strongly yᵢ pulls its own fitted value.
Pseudo-inverse (X⁺)	The Moore–Penrose generalised inverse; X⁺y gives the minimum-norm least-squares solution even when X is rank-deficient.
QR decomposition	Factorisation X = QR with orthonormal Q and upper-triangular R, used to solve least squares without forming XᵀX.
Singular value decomposition (SVD)	Factorisation X = UΣVᵀ; singular values in Σ reveal rank and conditioning.
Condition number	Ratio of largest to smallest singular value (of X) or eigenvalue (of XᵀX); measures sensitivity of the solution to small changes in data.
Gradient descent	An iterative method that updates parameters in the direction of the negative gradient of the loss.
Learning rate (η)	The step size multiplying the gradient in each gradient descent update.
Epoch	One full pass over the training data during iterative training.
Feature scaling / standardisation	Transforming features to comparable scales, commonly to mean 0 and standard deviation 1.
One-hot encoding	Representing a categorical feature with k levels as indicator columns that are 1 for the row's level and 0 otherwise.
Dummy variable trap	Perfect collinearity caused by including all k dummies of a categorical feature together with an intercept.
Reference (baseline) level	The category whose dummy is dropped; other dummies' coefficients are differences relative to it.
Interaction term	A product of features (for example x₁·x₂) that lets the effect of one feature depend on another.
Polynomial term	A power of a feature (for example x²) used to model curvature while staying linear in the parameters.
Gauss–Markov assumptions	Linearity in parameters, full rank, zero conditional mean of errors, and homoscedastic uncorrelated errors.
BLUE	Best Linear Unbiased Estimator: the linear unbiased estimator with the smallest variance; OLS under Gauss–Markov.
Exogeneity	The condition E[ε ∣ X] = 0; features carry no information about the error. Its failure (endogeneity) biases OLS.
Homoscedasticity	Constant error variance across observations.
Heteroscedasticity	Error variance that changes across observations; invalidates classical standard errors.
Breusch–Pagan test	A test for heteroscedasticity that regresses squared residuals on the features.
Robust (HC) standard errors	Standard errors computed with a sandwich formula that remains valid under heteroscedasticity.
Q-Q plot	A plot of sample residual quantiles against theoretical normal quantiles to judge normality.
Autocorrelation	Correlation between errors of different observations, typically adjacent in time.
Durbin–Watson statistic	A statistic between 0 and 4 for lag-1 residual autocorrelation; about 2 indicates none.
Multicollinearity	Strong linear relationships among features that inflate coefficient variances.
Variance Inflation Factor (VIF)	1/(1 − R²_j); the factor by which the variance of β̂_j is inflated by correlation with other features.
Influential point	An observation whose removal changes the fitted coefficients or predictions substantially.
Cook's distance	A per-row influence measure combining residual size and leverage.
Studentised residual	A residual divided by its estimated standard deviation σ̂√(1 − hᵢᵢ).
Standard error (SE)	The estimated standard deviation of a coefficient's sampling distribution.
t-statistic	A coefficient estimate divided by its standard error; used to test whether the coefficient is zero.
p-value	The probability, under the null hypothesis and model assumptions, of a test statistic at least as extreme as the one observed.
Confidence interval	A range computed by a procedure that contains the true parameter in a stated fraction of repeated samples.
Prediction interval	A range expected to contain a new individual observation at given feature values; wider than the confidence interval for the mean.
F-test	A test of whether a group of coefficients (all slopes, or a nested subset) are jointly zero.
Degrees of freedom (residual)	n − p − 1: the number of observations minus the number of estimated parameters.
R²	1 − SS_res/SS_tot; the fraction of variance in y explained by the model on the data evaluated.
Adjusted R²	R² penalised for the number of features: 1 − (1 − R²)(n − 1)/(n − p − 1).
RMSE	Root mean squared error, √(SS_res/n), in units of y.
MAE	Mean absolute error, (1/n)Σ∣yᵢ − ŷᵢ∣, in units of y.
Bias (of a model)	Systematic error from the model's inability to represent the true relationship.
Variance (of a model)	Variability of the model's predictions across different training samples.
Overfitting	Fitting noise in the training data so that performance on new data is worse than on training data.
Train/test split	Partitioning data into a part used for fitting and a held-out part used once for final evaluation.
k-fold cross-validation	Rotating each of k folds as validation data while training on the rest, then averaging scores.
Regularisation	Adding a penalty on coefficient size to the loss to reduce variance.
Ridge regression (L2)	Least squares plus λ‖β‖²; shrinks coefficients smoothly and always has a unique solution for λ > 0.
Lasso (L1)	Least squares plus a multiple of ‖β‖₁; shrinks and can set coefficients exactly to zero.
Elastic Net	A penalty mixing L1 and L2; sparse like lasso, stable with correlated features like ridge.
Soft-thresholding	The operation sign(z)·max(∣z∣ − λ, 0) that produces lasso solutions under orthonormal design.
Hyperparameter	A setting chosen before fitting (for example λ), usually tuned by cross-validation.
Data leakage	Use of information in training or validation that would not be available at prediction time.
Pipeline	A chained sequence of preprocessing steps and a final estimator that is fit and applied as one object.
Data drift	A change over time in the distribution of input features.
Concept drift	A change over time in the relationship between features and the target.
Training–serving skew	A mismatch between how features are computed during training and during serving.

---

Worked Example (a) — Multiple regression by hand with the normal equation

Setting (synthetic, for teaching): monthly cloud bill y (in units of $100) for 6 teams, explained by x₁ = number of VMs and x₂ = storage in TB.

Row	x₁ (VMs)	x₂ (TB)	y (bill, $100s)
1	1	2	5
2	2	1	7
3	3	3	12
4	4	3	13
5	5	5	19
6	6	4	19

n = 6, p = 2. Correlation between x₁ and x₂ is 0.8315 (correlated, but not perfectly).

Step 1 — Design matrix

      ⎡ 1  1  2 ⎤          ⎡  5 ⎤
      ⎢ 1  2  1 ⎥          ⎢  7 ⎥
  X = ⎢ 1  3  3 ⎥      y = ⎢ 12 ⎥
      ⎢ 1  4  3 ⎥          ⎢ 13 ⎥
      ⎢ 1  5  5 ⎥          ⎢ 19 ⎥
      ⎣ 1  6  4 ⎦          ⎣ 19 ⎦

Step 2 — XᵀX

Entries are sums: [n, Σx₁, Σx₂; Σx₁, Σx₁², Σx₁x₂; Σx₂, Σx₁x₂, Σx₂²].

Σx₁ = 21     Σx₂ = 18     Σx₁² = 91     Σx₂² = 64     Σx₁x₂ = 2+2+9+12+25+24 = 74

        ⎡  6  21  18 ⎤
XᵀX  =  ⎢ 21  91  74 ⎥
        ⎣ 18  74  64 ⎦

Step 3 — (XᵀX)⁻¹ via determinant and adjugate

det(XᵀX) = 324

               ⎡ 348   −12   −84 ⎤
adj(XᵀX)   =   ⎢ −12    60   −66 ⎥
               ⎣ −84   −66   105 ⎦

                               ⎡ 29/27    −1/27    −7/27  ⎤     ⎡  1.074074  −0.037037  −0.259259 ⎤
(XᵀX)⁻¹ = adj(XᵀX) / 324   =   ⎢ −1/27     5/27   −11/54  ⎥  ≈  ⎢ −0.037037   0.185185  −0.203704 ⎥
                               ⎣ −7/27   −11/54    35/108 ⎦     ⎣ −0.259259  −0.203704   0.324074 ⎦

Step 4 — Xᵀy

Σy = 75     Σx₁y = 5+14+36+52+95+114 = 316     Σx₂y = 10+7+36+39+95+76 = 263

Xᵀy = [ 75, 316, 263 ]ᵀ

Step 5 — β̂ = (XᵀX)⁻¹Xᵀy

adj(XᵀX)·Xᵀy = [ 216, 702, 459 ]ᵀ          (e.g. 348·75 − 12·316 − 84·263 = 216)

β̂ = [ 216, 702, 459 ]ᵀ / 324 = [ 2/3, 13/6, 17/12 ]ᵀ ≈ [ 0.6667, 2.1667, 1.4167 ]ᵀ

Fitted model: ŷ = 0.6667 + 2.1667·x₁ + 1.4167·x₂

Interpretation: holding storage constant, each additional VM is associated with an increase of about 2.1667 in the predicted bill ($216.67). Holding the number of VMs constant, each additional TB is associated with about 1.4167 ($141.67). The intercept (0.6667, i.e. $66.67) is the predicted bill at 0 VMs and 0 TB, which lies outside the data range and is a mathematical anchor, not a reliable estimate.

Step 6 — Predictions and residuals

Row	y	ŷ (exact)	ŷ	e = y − ŷ	leverage hᵢᵢ	Cook's D
1	5	17/3	5.6667	−0.6667	0.6296	1.1657
2	7	77/12	6.4167	+0.5833	0.6574	1.0891
3	12	137/12	11.4167	+0.5833	0.2130	0.0669
4	13	163/12	13.5833	−0.5833	0.2130	0.0669
5	19	223/12	18.5833	+0.4167	0.6574	0.5557
6	19	58/3	19.3333	−0.3333	0.6296	0.2914

Checks: Σeᵢ = 0 and Xᵀe = [0, 0, 0]ᵀ exactly (orthogonality from Section 3). Σhᵢᵢ = 3 = p + 1.

Step 7 — Sums of squares, R², adjusted R², RMSE, MAE

ȳ = 75/6 = 12.5
SS_res = Σe²        = 7/4   = 1.75
SS_tot = Σ(y − ȳ)²  = 343/2 = 171.5        (deviations −7.5, −5.5, −0.5, 0.5, 6.5, 6.5)
SS_reg = Σ(ŷ − ȳ)²  = 679/4 = 169.75       (check: 169.75 + 1.75 = 171.5)

R²          = 1 − 1.75/171.5            = 97/98   ≈ 0.98980
adjusted R² = 1 − (1 − 97/98)·(5/3)     = 289/294 ≈ 0.98299
RMSE        = √(1.75/6)                 ≈ 0.5401   ($54.01)
MAE         = (1/6)·Σ|e|                ≈ 0.5278

Step 8 — Inference (verified with statsmodels)

σ̂² = SS_res/(n − p − 1) = 1.75/3 = 7/12 ≈ 0.5833      σ̂ ≈ 0.7638
SE(β̂_j) = σ̂·√[(XᵀX)⁻¹]_jj

Coefficient	Estimate	SE	t	p-value (df = 3)	95% CI
β₀ (intercept)	0.6667	0.7915	0.842	0.4615	[−1.8524, 3.1857]
β₁ (VMs)	2.1667	0.3287	6.592	0.0071	[1.1207, 3.2126]
β₂ (TB)	1.4167	0.4348	3.258	0.0472	[0.0330, 2.8004]

t_{0.975, 3} = 3.1824. Overall F = (169.75/2)/(1.75/3) = 145.5 with (2, 3) df, p ≈ 0.00103.

Prediction at x₁ = 4, x₂ = 4: ŷ = 2/3 + 4·13/6 + 4·17/12 = 15.0. 95% CI for the mean: [13.5967, 16.4033]. 95% prediction interval for a new team: [12.1933, 17.8067].

Critical reading: with n = 6 and 3 parameters there are only 3 residual degrees of freedom, so p-values and intervals are fragile. Rows 1 and 2 have Cook's D above 1 because they sit at the edge of the feature space (leverage ≈ 0.63–0.66 versus an average of 3/6 = 0.5); with so few rows almost every point is influential. This is a hand-computation example, not a model to deploy.

Worked Example (b) — Multicollinearity and VIF

Setting (synthetic): throughput y (thousand requests/min) for 8 service configurations; x₁ = vCPUs, x₂ = memory in GB (provisioned at roughly 2 GB per vCPU, so nearly a copy of x₁), x₃ = number of cache nodes (varies independently).

Row	x₁	x₂	x₃	y
1	1	2.1	5	10.4
2	2	3.9	3	13.1
3	3	6.2	6	19.2
4	4	7.8	2	19.9
5	5	10.1	7	27.0
6	6	12.2	4	29.4
7	7	13.8	8	34.9
8	8	16.1	1	34.6

Step 1 — Correlations and auxiliary regressions

corr(x₁, x₂) = 0.99942, corr(x₁, x₃) = −0.0476, corr(x₂, x₃) = −0.0465.

For each feature, regress it on the other two (with intercept) and compute VIF_j = 1/(1 − R²_j):

Feature	R²_j	VIF_j	SE inflation √VIF
x₁ (vCPU)	0.998841	862.59	≈ 29.4
x₂ (memory)	0.998841	862.50	≈ 29.4
x₃ (cache)	0.003376	1.0034	≈ 1.00

statsmodels' variance_inflation_factor returns the same values (862.5922, 862.4977, 1.0034).

Step 2 — Symptoms in the OLS fit

Term	Coefficient	SE	p-value
Intercept	3.6735	0.2537	0.00013
x₁	1.1764	1.0763	0.3358
x₂	1.3096	0.5384	0.0718
x₃	0.6169	0.0367	0.00007

R² = 0.99963 and the overall F-test p-value is about 2.6 × 10⁻⁷, yet neither x₁ nor x₂ is individually significant at 0.05. This pattern (strong joint fit, weak individual t-tests) is a classic sign of collinearity. x₃, which is not collinear, has a small SE.

Condition numbers: statsmodels' Cond. No. on the unscaled design with intercept is 172.9; on standardised features cond(Z) = 58.79 and cond(ZᵀZ) = 3455.78. Eigenvalues of the feature correlation matrix are 0.00058, 0.9956, 2.0038: the near-zero eigenvalue identifies the near-dependence x₂ ≈ 2·x₁.

Step 3 — Instability under a tiny data change

Change one target value, row 4: y = 19.9 → 20.9 (one unit), and refit.

Term	Original fit	After change
x₁	1.1764	3.4254
x₂	1.3096	0.1764
x₃	0.6169	0.5593

The individual coefficients of x₁ and x₂ swing dramatically, while the combined effect along the shared direction, β₁ + 2·β₂, barely moves (3.7956 → 3.7782). The data can estimate the joint effect of "more vCPU with proportional memory" but cannot separate the two.

Remedy shown: drop x₂. The reduced model gives x₁ coefficient 3.7926 (SE 0.0517) with R² = 0.99909, and the same perturbation only moves it to 3.7778. Alternative remedy: ridge (Example c).

Worked Example (c) — Ridge vs OLS on small, collinear data

Same 8 rows as Example (b). Features are standardised (mean 0, population SD 1, as StandardScaler does), y is centered, and the intercept (= ȳ) is not penalised. Ridge coefficients are computed with the closed form w = (ZᵀZ + λI)⁻¹Zᵀ(y − ȳ) and cross-checked against sklearn.linear_model.Ridge(alpha=λ) (they match to floating-point precision).

Coefficients on the standardised scale (change in y per 1 SD of the feature)

λ	w₁ (vCPU)	w₂ (memory)	w₃ (cache)	‖w‖	train R²	cond(ZᵀZ + λI)
0 (OLS)	2.6955	5.9976	1.4135	6.7257	0.9996	3455.78
0.1	4.2451	4.3932	1.3953	6.2664	0.9995	154.16
1	4.0784	4.0953	1.2364	5.9105	0.9958	16.95
10	2.6610	2.6633	0.5587	3.8061	0.8451	2.60

Intercept for all rows: ȳ = 23.5625.

Same fits after the one-unit change to row 4 (from Example b)

λ	w₁	w₂	w₃
0 (OLS)	7.8487	0.8079	1.2814
0.1	4.4561	4.1456	1.2592
1	4.0854	4.0544	1.1138
10	2.6532	2.6506	0.4976

Leave-one-out cross-validated RMSE (scaler re-fit inside each fold)

Model	LOOCV RMSE
OLS	0.3551
Ridge λ = 0.1	0.2908
Ridge λ = 1	0.8768
Ridge λ = 10	4.7007

What the numbers show

1. Stability: under OLS the split between w₁ and w₂ changes from (2.70, 6.00) to (7.85, 0.81) after a one-unit change in one target. With λ = 0.1 it changes only from (4.25, 4.39) to (4.46, 4.15). Ridge spreads weight evenly across the two near-duplicate features.
2. Conditioning: adding λ to the diagonal drops the condition number from 3455.78 to 154.16 (λ = 0.1) and 16.95 (λ = 1).
3. Bias–variance trade-off: training R² always decreases as λ grows (ridge cannot beat OLS on training data). Out-of-sample error first improves (λ = 0.1 has the lowest LOOCV RMSE here, 0.2908 vs 0.3551 for OLS), then worsens as shrinkage becomes too strong (λ = 10 underfits). λ must be tuned by cross-validation.
4. No exact zeros: ridge shrinks all coefficients but none becomes exactly 0. Lasso would be needed for that (Lab 3 shows a lasso zero).

All numbers in Worked Examples (a), (b) and (c) were computed and cross-checked in Python 3.13 with numpy 2.5.3, scikit-learn 1.9.1, statsmodels 0.15.0 and sympy (exact fractions for Example a). Scripts: verify/verify_examples.py, verify/verify_extra_a.py, verify/verify_c_loo.py, verify/verify_misc.py.

Common setup (do once)

Prerequisites: Python 3.10 or newer, pip, about 300 MB free disk space for packages.

mkdir -p ~/lr-expert-labs && cd ~/lr-expert-labs
python3 -m venv .venv
source .venv/bin/activate            # Windows PowerShell: .venv\Scripts\Activate.ps1
pip install --upgrade pip
pip install numpy scipy pandas scikit-learn statsmodels matplotlib
python -c "import numpy, sklearn, statsmodels; print(numpy.__version__, sklearn.__version__, statsmodels.__version__)"

Using a virtual environment keeps these packages isolated from the system Python, and cleanup becomes deleting one folder.

Lab 1 — Normal equation from scratch, cross-checked with lstsq, QR, pinv and scikit-learn

Goal: implement β̂ = (XᵀX)⁻¹Xᵀy, confirm that six different routes give the same coefficients, verify the algebraic properties of the OLS solution (Xᵀe = 0, trace(H) = p + 1, H² = H), and measure why explicit inversion is avoided.

Prerequisites: common setup done; knowledge of Notes Sections 3–5.

Steps

1. Save the script below as lab1_normal_equation.py:

"""Lab 1: OLS via the normal equation from scratch, cross-checked with
np.linalg.lstsq, QR, pinv and scikit-learn LinearRegression."""
import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(42)
n, p = 200, 3
X_raw = rng.normal(size=(n, p))                     # 3 features
beta_true = np.array([4.0, 2.0, -3.0, 0.5])          # [intercept, b1, b2, b3]
X = np.column_stack([np.ones(n), X_raw])             # design matrix with intercept column
y = X @ beta_true + rng.normal(scale=1.0, size=n)    # noise sd = 1

# 1) Normal equation, two ways
XtX, Xty = X.T @ X, X.T @ y
beta_inv   = np.linalg.inv(XtX) @ Xty                # explicit inverse (teaching only)
beta_solve = np.linalg.solve(XtX, Xty)               # solve linear system (better)

# 2) Library / factorisation routes
beta_lstsq, _, rank, sv = np.linalg.lstsq(X, y, rcond=None)   # SVD-based
Q, R = np.linalg.qr(X)
beta_qr   = np.linalg.solve(R, Q.T @ y)              # QR route
beta_pinv = np.linalg.pinv(X) @ y                    # Moore-Penrose pseudo-inverse
sk = LinearRegression().fit(X_raw, y)                # sklearn adds intercept itself
beta_sk = np.r_[sk.intercept_, sk.coef_]

print("true beta      :", beta_true)
for name, b in [("inv", beta_inv), ("solve", beta_solve), ("lstsq", beta_lstsq),
                ("qr", beta_qr), ("pinv", beta_pinv), ("sklearn", beta_sk)]:
    print(f"{name:<15}: {np.round(b, 4)}   max|diff vs lstsq| = {np.max(np.abs(b - beta_lstsq)):.2e}")
print("rank(X) =", rank, " singular values:", np.round(sv, 3))

# 3) Properties of the OLS solution
y_hat = X @ beta_lstsq
e = y - y_hat
H = X @ np.linalg.solve(XtX, X.T)                    # hat matrix
print("sum of residuals          :", f"{e.sum():.2e}")
print("max |X^T e| (orthogonality):", f"{np.max(np.abs(X.T @ e)):.2e}")
print("trace(H) (= p+1)           :", round(np.trace(H), 6))
print("H idempotent (max|HH-H|)   :", f"{np.max(np.abs(H @ H - H)):.2e}")
r2 = 1 - (e @ e) / np.sum((y - y.mean())**2)
print("R^2 (manual) =", round(r2, 4), " R^2 (sklearn) =", round(sk.score(X_raw, y), 4))

# 4) Why not invert X^T X?  Condition number squares.
x1 = rng.normal(size=n)
X_bad = np.column_stack([np.ones(n), x1, x1 + 1e-6 * rng.normal(size=n)])
print("cond(X_bad)       = %.3e" % np.linalg.cond(X_bad))
print("cond(X_bad^T X_bad) = %.3e  (approx. cond(X)^2)" % np.linalg.cond(X_bad.T @ X_bad))

2. Run it:

python lab1_normal_equation.py

Expected output

true beta      : [ 4.   2.  -3.   0.5]
inv            : [ 3.9629  2.0368 -3.0076  0.4483]   max|diff vs lstsq| = 2.66e-15
solve          : [ 3.9629  2.0368 -3.0076  0.4483]   max|diff vs lstsq| = 2.22e-15
lstsq          : [ 3.9629  2.0368 -3.0076  0.4483]   max|diff vs lstsq| = 0.00e+00
qr             : [ 3.9629  2.0368 -3.0076  0.4483]   max|diff vs lstsq| = 1.33e-15
pinv           : [ 3.9629  2.0368 -3.0076  0.4483]   max|diff vs lstsq| = 6.66e-15
sklearn        : [ 3.9629  2.0368 -3.0076  0.4483]   max|diff vs lstsq| = 1.78e-15
rank(X) = 4  singular values: [14.665 14.421 13.663 12.599]
sum of residuals          : -3.60e-13
max |X^T e| (orthogonality): 3.64e-13
trace(H) (= p+1)           : 4.0
H idempotent (max|HH-H|)   : 4.16e-17
R^2 (manual) = 0.921  R^2 (sklearn) = 0.921
cond(X_bad)       = 2.089e+06
cond(X_bad^T X_bad) = 4.358e+12  (approx. cond(X)^2)

Verification

• All six coefficient vectors agree to about 1e-14 or better (max|diff vs lstsq|). The estimates (3.96, 2.04, −3.01, 0.45) are close to, but not equal to, the true values (4, 2, −3, 0.5) because of the noise term.
• sum of residuals and max |X^T e| are zero up to round-off: the normal equations hold.
• trace(H) = 4.0 equals p + 1 (3 features + intercept).
• cond(X_bad^T X_bad) is roughly the square of cond(X_bad): (2.089e6)² ≈ 4.36e12. With about 16 significant digits in double precision, forming XᵀX here leaves only around 4 reliable digits, which is why QR/SVD are used.

Extension (optional)

Append these lines to the script to see an exactly singular design (the last two columns identical):

X_sing = np.column_stack([np.ones(n), x1, x1])
try:
    print("inv:", np.linalg.inv(X_sing.T @ X_sing) @ X_sing.T @ y)
except np.linalg.LinAlgError as err:
    print("inv failed:", err)
b, _, r, _ = np.linalg.lstsq(X_sing, y, rcond=None)
print("lstsq:", np.round(b, 4), "rank =", r)

Reference output:

inv failed: Singular matrix
lstsq: [4.1809 0.0067 0.0067] rank = 2

(y was generated without this new x1, so its coefficient is near 0.) In the reference run np.linalg.inv raised LinAlgError: Singular matrix (on other inputs it can instead return meaningless huge values without an error), while np.linalg.lstsq returned the minimum-norm solution, split the shared coefficient equally between the two identical columns, and reported rank = 2.

Cleanup

rm lab1_normal_equation.py

Lab 2 — Batch gradient descent with feature scaling, loss curve, match to OLS

Goal: implement batch gradient descent for MSE, show that unscaled features force a tiny learning rate (or diverge), show that standardised features converge fast, map the scaled coefficients back to original units, and match the closed-form OLS solution.

Prerequisites: common setup done; Lab 1; Notes Section 6.

Steps

1. Save the script below as lab2_gradient_descent.py:

"""Lab 2: Batch gradient descent for multiple linear regression, with and
without feature scaling, compared with the closed-form OLS solution."""
import numpy as np
import matplotlib
matplotlib.use("Agg")                     # no display needed; writes PNG
import matplotlib.pyplot as plt

rng = np.random.default_rng(7)
n = 500
size_sqft = rng.uniform(500, 3500, n)    # large scale feature
rooms     = rng.integers(1, 6, n).astype(float)
age_years = rng.uniform(0, 40, n)
X_raw = np.column_stack([size_sqft, rooms, age_years])
y = 50 + 0.12 * size_sqft + 8.0 * rooms - 1.5 * age_years + rng.normal(0, 15, n)

def mse(X, y, w):                         # X includes the intercept column
    r = X @ w - y
    return (r @ r) / len(y)

def batch_gd(X, y, lr, epochs):
    np.seterr(over="ignore", invalid="ignore")   # divergence is shown on purpose
    w = np.zeros(X.shape[1])
    hist = []
    for _ in range(epochs):
        grad = (2 / len(y)) * X.T @ (X @ w - y)   # gradient of MSE
        w -= lr * grad
        hist.append(mse(X, y, w))
        if not np.isfinite(hist[-1]):
            break
    return w, np.array(hist)

# Closed-form reference on the original scale
X1 = np.column_stack([np.ones(n), X_raw])
w_ols = np.linalg.lstsq(X1, y, rcond=None)[0]
print("OLS (lstsq) original scale :", np.round(w_ols, 4))

# A) Unscaled features
for lr in [1e-7, 1e-6]:
    w_u, h_u = batch_gd(X1, y, lr=lr, epochs=2000)
    print(f"unscaled lr={lr:g}: final MSE={h_u[-1]:.4g}  epochs run={len(h_u)}")

# B) Standardised features (fit mean/std on the data used for training)
mu, sd = X_raw.mean(0), X_raw.std(0)
Xs = np.column_stack([np.ones(n), (X_raw - mu) / sd])
w_s, h_s = batch_gd(Xs, y, lr=0.1, epochs=500)
# Map back to the original scale: b_j = w_j / sd_j ; b0 = w0 - sum(w_j * mu_j / sd_j)
b = w_s[1:] / sd
b0 = w_s[0] - np.sum(w_s[1:] * mu / sd)
w_back = np.r_[b0, b]
print("GD (scaled) mapped back     :", np.round(w_back, 4))
print("max |GD - OLS|              : %.2e" % np.max(np.abs(w_back - w_ols)))
print("final MSE (GD scaled)       : %.4f" % h_s[-1])
print("MSE of OLS                  : %.4f" % mse(X1, y, w_ols))
print("cond(X^T X), unscaled design : %.3e" % np.linalg.cond(X1.T @ X1))
print("cond(X^T X), scaled design   : %.3f" % np.linalg.cond(Xs.T @ Xs))
# GD on MSE is stable only if lr < 1 / lambda_max((1/n) X^T X)
for tag, A in [("unscaled", X1), ("scaled", Xs)]:
    lmax = np.linalg.eigvalsh(A.T @ A / n).max()
    print(f"largest stable lr, {tag:<8}: {1/lmax:.3e}")

# Too large a learning rate on scaled data
_, h_big = batch_gd(Xs, y, lr=1.1, epochs=50)
print("scaled lr=1.1: MSE after %d epochs = %.3e (diverges)" % (len(h_big), h_big[-1]))

plt.figure(figsize=(7, 4))
plt.plot(h_s, label="scaled, lr=0.1")
plt.yscale("log"); plt.xlabel("epoch"); plt.ylabel("training MSE (log scale)")
plt.axhline(mse(X1, y, w_ols), ls="--", c="k", label="OLS minimum")
plt.legend(); plt.title("Batch gradient descent loss curve"); plt.tight_layout()
plt.savefig("lab2_loss_curve.png", dpi=120)
print("saved lab2_loss_curve.png")

2. Run it:

python lab2_gradient_descent.py

Expected output

OLS (lstsq) original scale : [54.4472  0.1189  7.4043 -1.513 ]
unscaled lr=1e-07: final MSE=1017  epochs run=2000
unscaled lr=1e-06: final MSE=inf  epochs run=160
GD (scaled) mapped back     : [54.4472  0.1189  7.4043 -1.513 ]
max |GD - OLS|              : 7.82e-14
final MSE (GD scaled)       : 245.0873
MSE of OLS                  : 245.0873
cond(X^T X), unscaled design : 6.667e+07
cond(X^T X), scaled design   : 1.056
largest stable lr, unscaled: 2.059e-07
largest stable lr, scaled  : 9.745e-01
scaled lr=1.1: MSE after 50 epochs = 3.573e+13 (diverges)
saved lab2_loss_curve.png

A file lab2_loss_curve.png is written in the current folder. It shows the training MSE (log scale) dropping from about 6 × 10⁴ to the OLS minimum (dashed line, MSE = 245.09) within roughly the first 20–30 epochs and then staying flat.

Verification

• GD (scaled) mapped back equals OLS (lstsq) and max |GD - OLS| is below 1e-10.
• The computed stability limit explains the unscaled runs: lr = 1e-7 is below the limit 2.06e-7 so it is stable, but after 2000 epochs MSE is still about 1017 (far from 245) because convergence is slow when cond(XᵀX) ≈ 6.7e7. lr = 1e-6 is above the limit and diverges to inf.
• With scaling, cond(XᵀX) ≈ 1.06, the limit is about 0.97, and lr = 0.1 converges; lr = 1.1 exceeds the limit and diverges.
• Open lab2_loss_curve.png and confirm the curve reaches the dashed line.

Cleanup

rm lab2_gradient_descent.py lab2_loss_curve.png

Lab 3 — End-to-end multiple regression on the scikit-learn diabetes dataset

Dataset: sklearn.datasets.load_diabetes: 442 patients, 10 baseline features (age, sex, body mass index bmi, average blood pressure bp, six blood serum measurements s1–s6), target = a quantitative measure of disease progression one year after baseline. The data file ships inside the scikit-learn package, so it loads offline. The lab uses scaled=False to get features in original units (by default scikit-learn returns features already mean-centered and scaled, which would hide the role of our own scaler). sex is coded 1/2 in the raw data and is used as a numeric binary feature.

Goal: a leakage-safe Pipeline (scaling + model), train/test split, 5-fold CV, statsmodels OLS inference, VIF, residual diagnostics (Breusch–Pagan, Q-Q, Durbin–Watson, influence), heteroscedasticity-robust standard errors, and a RidgeCV / LassoCV comparison on the same split.

Prerequisites: common setup done; Labs 1–2; Notes Sections 10–16.

Steps

1. Save the script below as lab3_diabetes_end_to_end.py:

"""Lab 3: End-to-end multiple linear regression on the scikit-learn
diabetes dataset (ships inside scikit-learn; no download needed)."""
import numpy as np, pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.stattools import durbin_watson, jarque_bera
from scipy import stats
from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split, KFold, cross_val_score
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error

pd.set_option("display.width", 120)
# ---- 1. Load raw (unscaled) features so that our own scaler matters --------
data = load_diabetes(as_frame=True, scaled=False)
X, y = data.data, data.target
print("Step 1 | shape:", X.shape, "| target mean = %.2f, sd = %.2f" % (y.mean(), y.std()))

# ---- 2. Hold out a test set BEFORE any fitting (prevents leakage) ----------
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)
print("Step 2 | train:", X_tr.shape, " test:", X_te.shape)

# ---- 3. Pipeline + 5-fold CV on the training set only ----------------------
kf = KFold(n_splits=5, shuffle=True, random_state=42)
ols_pipe = Pipeline([("scale", StandardScaler()), ("model", LinearRegression())])
cv_r2 = cross_val_score(ols_pipe, X_tr, y_tr, cv=kf, scoring="r2")
cv_rmse = -cross_val_score(ols_pipe, X_tr, y_tr, cv=kf, scoring="neg_root_mean_squared_error")
print("Step 3 | OLS 5-fold CV R2  : mean %.4f  sd %.4f  folds %s" % (cv_r2.mean(), cv_r2.std(), np.round(cv_r2, 3)))
print("       | OLS 5-fold CV RMSE: mean %.2f" % cv_rmse.mean())

def evaluate(name, pipe):
    pipe.fit(X_tr, y_tr)
    pred = pipe.predict(X_te)
    return dict(model=name,
                train_R2=r2_score(y_tr, pipe.predict(X_tr)),
                test_R2=r2_score(y_te, pred),
                test_RMSE=np.sqrt(mean_squared_error(y_te, pred)),
                test_MAE=mean_absolute_error(y_te, pred))

rows = [evaluate("OLS", ols_pipe)]
print("Step 3 | OLS test R2 = %.4f, RMSE = %.2f" % (rows[0]["test_R2"], rows[0]["test_RMSE"]))

# ---- 4. statsmodels OLS for inference (training data, original units) ------
Xc_tr = sm.add_constant(X_tr)
sm_fit = sm.OLS(y_tr, Xc_tr).fit()
print("\nStep 4 | statsmodels OLS summary (training set)")
print(sm_fit.summary())

# ---- 5. VIF (training features, with constant in the design) ---------------
vif = pd.Series([variance_inflation_factor(Xc_tr.values, i) for i in range(1, Xc_tr.shape[1])],
                index=X_tr.columns, name="VIF").round(2)
print("\nStep 5 | VIF per feature\n", vif.sort_values(ascending=False).to_string())

# ---- 6. Residual diagnostics ----------------------------------------------
resid, fitted = sm_fit.resid, sm_fit.fittedvalues
bp_lm, bp_p, _, _ = het_breuschpagan(resid, Xc_tr)
jb, jb_p, skew, kurt = jarque_bera(resid)
sh_w, sh_p = stats.shapiro(resid)
dw = durbin_watson(resid)
infl = sm_fit.get_influence()
cooks = infl.cooks_distance[0]
lev = infl.hat_matrix_diag
n_tr, k = Xc_tr.shape
print("\nStep 6 | Breusch-Pagan LM = %.3f, p = %.4f" % (bp_lm, bp_p))
print("       | Jarque-Bera = %.3f, p = %.4f (skew %.3f, kurtosis %.3f)" % (jb, jb_p, skew, kurt))
print("       | Shapiro-Wilk W = %.4f, p = %.4f" % (sh_w, sh_p))
print("       | Durbin-Watson = %.3f" % dw)
print("       | points with Cook's D > 4/n (%.4f): %d ; max Cook's D = %.4f" % (4 / n_tr, (cooks > 4 / n_tr).sum(), cooks.max()))
print("       | points with leverage > 2k/n (%.4f): %d" % (2 * k / n_tr, (lev > 2 * k / n_tr).sum()))
print("       | statsmodels condition number (unscaled design) = %.1f" % sm_fit.condition_number)

# Heteroscedasticity-robust (HC3) standard errors: same coefficients, different SEs
hc3 = sm.OLS(y_tr, Xc_tr).fit(cov_type="HC3")
cmp = pd.DataFrame({"coef_OLS": sm_fit.params, "coef_HC3": hc3.params,
                    "SE_classic": sm_fit.bse, "SE_HC3": hc3.bse,
                    "p_classic": sm_fit.pvalues, "p_HC3": hc3.pvalues}).round(4)
print("       | classic vs HC3 robust inference (coefficients identical, SEs differ)")
print(cmp.to_string())

fig, ax = plt.subplots(1, 2, figsize=(11, 4.2))
ax[0].scatter(fitted, resid, s=12, alpha=.6); ax[0].axhline(0, c="k", lw=1)
ax[0].set_xlabel("fitted value"); ax[0].set_ylabel("residual"); ax[0].set_title("Residuals vs fitted")
sm.qqplot(resid, line="45", fit=True, ax=ax[1]); ax[1].set_title("Normal Q-Q of residuals")
plt.tight_layout(); plt.savefig("lab3_residual_diagnostics.png", dpi=120)
print("       | saved lab3_residual_diagnostics.png")

# ---- 7. Regularised models: RidgeCV and LassoCV inside pipelines -----------
ridge_pipe = Pipeline([("scale", StandardScaler()),
                       ("model", RidgeCV(alphas=np.logspace(-3, 3, 61)))])
lasso_pipe = Pipeline([("scale", StandardScaler()),
                       ("model", LassoCV(cv=kf, random_state=42, max_iter=50000))])
rows.append(evaluate("RidgeCV", ridge_pipe))
rows.append(evaluate("LassoCV", lasso_pipe))
print("\nStep 7 | RidgeCV alpha = %.4f" % ridge_pipe.named_steps["model"].alpha_)
print("       | LassoCV alpha = %.4f" % lasso_pipe.named_steps["model"].alpha_)

coefs = pd.DataFrame({
    "OLS": ols_pipe.named_steps["model"].coef_,
    "Ridge": ridge_pipe.named_steps["model"].coef_,
    "Lasso": lasso_pipe.named_steps["model"].coef_}, index=X.columns).round(2)
print("       | coefficients on the standardised scale (target units per 1 SD of feature)")
print(coefs.to_string())
print("       | Lasso zero coefficients:", list(coefs.index[coefs["Lasso"] == 0]))

print("\nStep 8 | Test-set comparison (same split, random_state=42)")
print(pd.DataFrame(rows).set_index("model").round(4).to_string())

2. Run it:

python lab3_diabetes_end_to_end.py

Expected output

Step 1 | shape: (442, 10) | target mean = 152.13, sd = 77.09
Step 2 | train: (353, 10)  test: (89, 10)
Step 3 | OLS 5-fold CV R2  : mean 0.4804  sd 0.0409  folds [0.47  0.537 0.411 0.491 0.493]
       | OLS 5-fold CV RMSE: mean 55.39
Step 3 | OLS test R2 = 0.4526, RMSE = 53.85

Step 4 | statsmodels OLS summary (training set)
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 target   R-squared:                       0.528
Model:                            OLS   Adj. R-squared:                  0.514
Method:                 Least Squares   F-statistic:                     38.25
Date:                Fri, 02 Oct 2026   Prob (F-statistic):           5.41e-50
Time:                        04:48:13   Log-Likelihood:                -1906.1
No. Observations:                 353   AIC:                             3834.
Df Residuals:                     342   BIC:                             3877.
Df Model:                          10                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const       -341.3782     74.073     -4.609      0.000    -487.074    -195.683
age            0.1377      0.251      0.549      0.583      -0.356       0.631
sex          -23.0645      6.536     -3.529      0.000     -35.921     -10.208
bmi            5.8464      0.829      7.049      0.000       4.215       7.478
bp             1.1971      0.246      4.873      0.000       0.714       1.680
s1            -1.2817      0.621     -2.065      0.040      -2.503      -0.061
s2             0.8112      0.570      1.423      0.156      -0.310       1.933
s3             0.6017      0.858      0.701      0.484      -1.086       2.289
s4            10.1595      6.841      1.485      0.138      -3.297      23.616
s5            67.1090     17.542      3.826      0.000      32.606     101.612
s6             0.2016      0.304      0.663      0.508      -0.397       0.800
==============================================================================
Omnibus:                        1.457   Durbin-Watson:                   1.794
Prob(Omnibus):                  0.483   Jarque-Bera (JB):                1.412
Skew:                           0.064   Prob(JB):                        0.494
Kurtosis:                       2.718   Cond. No.                     7.07e+03
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 7.07e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

Step 5 | VIF per feature
 s1     55.25
s2     35.76
s3     14.29
s5     10.07
s4      9.33
bmi     1.57
s6      1.50
bp      1.42
sex     1.27
age     1.22

Step 6 | Breusch-Pagan LM = 19.484, p = 0.0345
       | Jarque-Bera = 1.412, p = 0.4937 (skew 0.064, kurtosis 2.718)
       | Shapiro-Wilk W = 0.9965, p = 0.6398
       | Durbin-Watson = 1.794
       | points with Cook's D > 4/n (0.0113): 18 ; max Cook's D = 0.0277
       | points with leverage > 2k/n (0.0623): 22
       | statsmodels condition number (unscaled design) = 7065.5
       | classic vs HC3 robust inference (coefficients identical, SEs differ)
       coef_OLS  coef_HC3  SE_classic   SE_HC3  p_classic   p_HC3
const -341.3782 -341.3782     74.0728  74.2562     0.0000  0.0000
age      0.1377    0.1377      0.2508   0.2485     0.5834  0.5795
sex    -23.0645  -23.0645      6.5362   6.3047     0.0005  0.0003
bmi      5.8464    5.8464      0.8294   0.8732     0.0000  0.0000
bp       1.1971    1.1971      0.2457   0.2556     0.0000  0.0000
s1      -1.2817   -1.2817      0.6207   0.5944     0.0397  0.0311
s2       0.8112    0.8112      0.5701   0.5336     0.1557  0.1285
s3       0.6017    0.6017      0.8579   0.8345     0.4836  0.4709
s4      10.1595   10.1595      6.8415   6.9594     0.1385  0.1443
s5      67.1090   67.1090     17.5418  16.9635     0.0002  0.0001
s6       0.2016    0.2016      0.3042   0.2916     0.5079  0.4893
       | saved lab3_residual_diagnostics.png

Step 7 | RidgeCV alpha = 1.2589
       | LassoCV alpha = 0.5420
       | coefficients on the standardised scale (target units per 1 SD of feature)
       OLS  Ridge  Lasso
age   1.75   1.82   1.29
sex -11.51 -11.43 -10.30
bmi  25.61  25.75  26.23
bp   16.83  16.71  16.13
s1  -44.45 -32.84 -11.37
s2   24.64  15.64   0.00
s3    7.68   2.57  -6.65
s4   13.14  11.51   7.12
s5   35.16  30.67  23.05
s6    2.35   2.48   2.31
       | Lasso zero coefficients: ['s2']

Step 8 | Test-set comparison (same split, random_state=42)
         train_R2  test_R2  test_RMSE  test_MAE
model                                          
OLS        0.5279   0.4526    53.8534   42.7941
RidgeCV    0.5275   0.4544    53.7650   42.8152
LassoCV    0.5244   0.4613    53.4261   42.8364

A file lab3_residual_diagnostics.png is written: left panel residuals vs fitted values (a band around 0 roughly −150 to +160 with no strong curve), right panel a normal Q-Q plot with points close to the 45° line.

Reading the results

• Generalisation: OLS 5-fold CV R² on the training set is 0.4804 (fold range 0.411–0.537); test R² is 0.4526 with RMSE 53.85. Train R² (0.5279) is only modestly above test R², so the model is not heavily overfitting; it is limited by the information in 10 features (bias), not by variance.
• Inference (Step 4): bmi, bp, s5 and sex have p < 0.001; s1 has p ≈ 0.040. age, s2, s3, s4, s6 are not individually significant. Adjusted R² (0.514) is below R² (0.528) because of the 10-feature penalty. The statsmodels condition-number warning (7.07e+03) is computed on unscaled columns with very different units; use VIF to judge collinearity.
• VIF (Step 5): s1 (55.25), s2 (35.76), s3 (14.29) and s5 (10.07) exceed the common VIF > 10 rule of thumb; s4 (9.33) exceeds 5. These serum measures are algebraically related (for example s4 is total cholesterol divided by HDL), so their individual coefficients and signs are unstable; this is why s1 and s2 have large opposite-signed standardised coefficients.
• Heteroscedasticity (Step 6): Breusch–Pagan p ≈ 0.0345 gives some evidence of non-constant variance at the 5% level. Coefficients are unchanged under HC3, but standard errors shift (for example bmi 0.829 → 0.873, s5 17.54 → 16.96). Conclusions about which features are significant at 0.05 do not change here.
• Normality: Jarque–Bera p ≈ 0.49 and Shapiro–Wilk p ≈ 0.64 show no evidence against normal residuals; the Q-Q plot agrees.
• Autocorrelation: Durbin–Watson = 1.794. Rows are patients, not a time series, so DW mainly checks that there is no ordering artefact; a value near 2 is consistent with none.
• Influence: 18 rows exceed the Cook's D > 4/n rule of thumb, but the maximum is only 0.0277 (far below 1), so no single patient dominates the fit.
• Regularisation (Steps 7–8): RidgeCV chose alpha ≈ 1.26 and LassoCV alpha ≈ 0.54. Lasso set s2 exactly to 0 and shrank s1 from −44.45 to −11.37; ridge shrank the collinear serum coefficients but kept all non-zero. Test RMSE: OLS 53.85, Ridge 53.77, Lasso 53.43. The differences (under 0.5 on a target with SD 77) are small relative to the variability across CV folds, so on this split regularisation mainly buys stability and a simpler model, not a large accuracy gain. Do not over-interpret a single train/test split.

Verification

1. The test-set table shows three rows; Lasso test R² (0.4613) ≥ Ridge (0.4544) ≥ OLS (0.4526) on this split.
2. coef_OLS and coef_HC3 columns are identical; only standard errors and p-values differ.
3. Re-run with random_state=0 in train_test_split to see split variance. In the reference run the test R² dropped to 0.3322 (OLS), 0.3348 (Ridge) and 0.3354 (Lasso), from about 0.45–0.46 with random_state=42. The ranking of models is the same, but a single split moves R² by more than 0.1, which is why CV mean and spread matter more than one test number.
4. Optional leakage check: move StandardScaler outside the pipeline, fit it on all of X, then run CV. For plain OLS the scores do not change (OLS is scale-equivariant), but for RidgeCV/LassoCV the penalty would be tuned on statistics that include test rows. Keep scaling inside the pipeline.

Cleanup

rm lab3_diabetes_end_to_end.py lab3_residual_diagnostics.png

Full cleanup (all labs)

deactivate                      # leave the virtual environment
cd ~ && rm -rf ~/lr-expert-labs  # removes scripts, plots and the .venv folder

Nothing was created outside ~/lr-expert-labs; no cloud resources were used.

Interactive quiz — 10 advanced MCQs

Original teaching items (not from any real exam). Medium = 50 XP, Hard = 80 XP, +10 XP per consecutive correct answer after the first in a streak.

1. Cheat sheet

1.1 Model and estimator

Item	Formula
Model	y = Xβ + ε, X is n × (p+1) with a leading column of ones
OLS objective	minimise ‖y − Xβ‖²
Normal equations	XᵀXβ̂ = Xᵀy
OLS solution (full rank)	β̂ = (XᵀX)⁻¹Xᵀy — compute with QR/SVD (np.linalg.lstsq), not an explicit inverse
Hat matrix	H = X(XᵀX)⁻¹Xᵀ, ŷ = Hy, e = (I − H)y, trace(H) = p + 1, H² = H
Orthogonality	Xᵀe = 0; with intercept Σeᵢ = 0
Conditioning	cond(XᵀX) = cond(X)²
GD update (MSE)	β ← β − η·(2/n)·Xᵀ(Xβ − y); stable if η < 1/λ_max(XᵀX/n)
Back-transform from standardised	β_j = w_j/s_j, β₀ = w₀ − Σ w_j μ_j/s_j
Interaction effect	∂E[y]/∂x₁ = β₁ + β₃x₂ for β₃x₁x₂

1.2 Inference

Item	Formula
Var(β̂)	σ²(XᵀX)⁻¹ (homoscedastic, uncorrelated errors)
σ̂²	SS_res/(n − p − 1)
SE(β̂_j)	σ̂·√[(XᵀX)⁻¹]_jj
t-statistic	β̂_j / SE(β̂_j), df = n − p − 1
CI for β_j	β̂_j ± t_{1−α/2, n−p−1}·SE(β̂_j)
CI for mean at x₀	ŷ₀ ± t·σ̂·√(x₀(XᵀX)⁻¹x₀ᵀ)
Prediction interval at x₀	ŷ₀ ± t·σ̂·√(1 + x₀(XᵀX)⁻¹x₀ᵀ)
Overall F	(SS_reg/p) / (SS_res/(n − p − 1))
Nested F (q restrictions)	[(SS_res,R − SS_res,F)/q] / [SS_res,F/(n − p − 1)]

1.3 Metrics

Metric	Formula	Note
R²	1 − SS_res/SS_tot	Never decreases in-sample when features are added (OLS with intercept); can be negative on test data
Adjusted R²	1 − (1 − R²)(n − 1)/(n − p − 1)	Can decrease when a feature is added (rises only if that feature's |t| > 1)
RMSE	√(SS_res/n)	Units of y; penalises large errors
Residual standard error	√(SS_res/(n − p − 1))	What statsmodels/R report as σ̂
MAE	(1/n)Σ∣eᵢ∣	Units of y; less sensitive to outliers
LOOCV (OLS shortcut)	(1/n)Σ(eᵢ/(1 − hᵢᵢ))²	No refitting needed

1.4 Regularisation

Method	Objective (scikit-learn scaling)	Key property
Ridge	‖y − Xw‖² + α‖w‖²	Closed form (XᵀX + αI)⁻¹Xᵀy on centered data; always unique for α > 0; no exact zeros
Lasso	(1/(2n))‖y − Xw‖² + α‖w‖₁	Exact zeros via soft-thresholding; unstable choice among correlated features
Elastic Net	(1/(2n))‖y − Xw‖² + αρ‖w‖₁ + (α(1 − ρ)/2)‖w‖²	Sparse + grouping of correlated features
Orthonormal-design intuition	ridge z/(1 + λ); lasso sign(z)·max(∣z∣ − λ, 0)	(for ½‖·‖² loss scaling)

Always: standardise features inside a Pipeline; do not penalise the intercept; tune α by CV on a log grid.

1.5 Diagnostics and rules of thumb

These thresholds are conventions used for screening, not hard rules. Always look at plots and at the practical consequences.

Check	Tool	Common rule of thumb
Non-linearity	Residuals vs fitted / vs each feature	Visible curve → add terms or transform
Heteroscedasticity	Funnel in residual plot; Breusch–Pagan (het_breuschpagan), White	Small p-value (for example < 0.05) suggests non-constant variance → robust (HC3) SEs, transform y, or WLS
Normality	Q-Q plot; Jarque–Bera; Shapiro–Wilk	Matters mainly for small-n inference and prediction intervals; large n makes tests flag trivial deviations
Autocorrelation	Durbin–Watson (time-ordered data); Breusch–Godfrey	About 2 = none; roughly < 1.5 or > 2.5 worth investigating; use tabulated bounds for a formal test
Multicollinearity	VIF = 1/(1 − R²_j)	VIF > 5 deserves attention; VIF > 10 often treated as serious
Conditioning	Condition number of standardised X	> 30 often cited (Belsley–Kuh–Welsch) as moderate/strong collinearity; statsmodels warns at > 1000 on the unscaled design
Leverage	hᵢᵢ	> 2(p+1)/n (or 3(p+1)/n) flags unusual feature values
Influence	Cook's distance	> 1 large; > 4/n worth a look
Outliers in y	Studentised residuals	|r| > 2 (or 3) worth a look

2. Common mistakes

1. Choosing models by training R². R² never decreases as features are added. Compare models by cross-validated error (or at least adjusted R²/AIC/BIC for in-sample comparisons).
2. Forgetting the intercept in statsmodels. sm.OLS(y, X) does not add a constant; without sm.add_constant(X) the model is forced through the origin, and the reported R² is uncentered and not comparable.
3. Falling into the dummy variable trap. Encoding all k levels plus an intercept makes XᵀX singular. Use k − 1 dummies (drop="first") for OLS with inference.
4. Misreading coefficients. A coefficient is a partial association holding the other included features constant; it changes when features are added or removed, is causal only under extra assumptions, and carries units, so raw magnitudes are not comparable across features (standardise first, and remember collinearity makes individual magnitudes unstable).
5. Thinking heteroscedasticity biases coefficients. It does not (if exogeneity holds); it breaks the classical standard errors. Fix inference with robust SEs; do not drop features because of it.
6. Treating VIF cut-offs as laws. VIF > 5 or > 10 are rules of thumb. High VIF matters when you need individual coefficients; it matters much less for pure prediction if the collinearity pattern is stable.
7. Data leakage through preprocessing. Fitting scalers, imputers, encoders or feature selection on all data before splitting or CV. Put them inside a Pipeline.
8. Regularising unscaled features. Ridge and lasso penalties depend on feature units; without standardisation, features with large numeric ranges are penalised less.
9. Random splits on time series or grouped data. Use TimeSeriesSplit or GroupKFold; otherwise the validation score is optimistic. Relatedly, p-values reported after stepwise or lasso selection are over-optimistic.
10. Extrapolating. Linear models will predict outside the training range without warning; validate input ranges at serving time.

3. Interview questions with model answers

Q1. Derive the OLS estimator in matrix form. When does it fail to exist?
Minimise S(β) = ‖y − Xβ‖² = yᵀy − 2βᵀXᵀy + βᵀXᵀXβ. The gradient −2Xᵀy + 2XᵀXβ is zero at XᵀXβ̂ = Xᵀy; the Hessian 2XᵀX is positive semi-definite, so this is a minimum. If X has full column rank, XᵀX is invertible and β̂ = (XᵀX)⁻¹Xᵀy is unique. It fails to be unique when columns are linearly dependent (dummy trap, duplicated features in different units) or when p + 1 > n; then use the pseudo-inverse (minimum-norm solution), drop redundant columns, or use ridge.

Q2. Why do libraries not compute (XᵀX)⁻¹ directly?
Forming XᵀX squares the condition number (cond(XᵀX) = cond(X)²), so precision can be lost before the inverse is even computed, and explicit inversion is slower and less stable than solving. QR (X = QR, solve Rβ = Qᵀy) works at the conditioning of X; SVD additionally handles rank deficiency and gives the minimum-norm solution. np.linalg.lstsq uses an SVD-based LAPACK routine and scikit-learn's LinearRegression calls SciPy's lstsq.

Q3. Explain the hat matrix and leverage.
H = X(XᵀX)⁻¹Xᵀ orthogonally projects y onto the column space of X: ŷ = Hy. It is symmetric and idempotent with trace p + 1. The diagonal hᵢᵢ is the leverage of row i: how far its feature values are from the feature centroid and how strongly yᵢ pulls ŷᵢ. Average leverage is (p+1)/n; residual variance is σ²(1 − hᵢᵢ), so high-leverage points have smaller raw residuals, which is why studentised residuals are used.

Q4. State the Gauss–Markov assumptions. Which are needed for unbiasedness, and where does normality come in?
Linearity in parameters, full column rank, E[ε | X] = 0, and Var(ε | X) = σ²I. The first three give unbiasedness; adding the fourth makes OLS BLUE. Normality is not part of the theorem; it makes t and F statistics exactly t- and F-distributed in finite samples. With large n, inference is approximately valid without it.

Q5. What does heteroscedasticity do to OLS, and how do you handle it?
Coefficients stay unbiased and consistent, but OLS is no longer efficient and the classical σ²(XᵀX)⁻¹ variance is wrong, so t-tests and CIs are unreliable. Detect with residual-vs-fitted plots and Breusch–Pagan or White tests. Handle with heteroscedasticity-consistent (HC0–HC3) standard errors, a variance-stabilising transform of y (for example log), or weighted least squares if the variance structure is known.

Q6. What is multicollinearity, how do you measure it, and when does it matter?
Strong linear dependence among features. It inflates Var(β̂_j) by VIF_j = 1/(1 − R²_j), making coefficients unstable and individual t-tests weak while overall fit and the F-test can stay strong. Measure with VIF (rules of thumb > 5, > 10) and the condition number of standardised X (> 30 often cited). It matters for interpreting individual coefficients; for prediction it matters less if the same collinearity pattern holds in future data. Remedies: drop/combine features, center before forming powers or interactions, ridge, PCA/PLS.

Q7. Ridge vs lasso vs Elastic Net — how do you choose?
Ridge (L2) shrinks all coefficients smoothly, handles collinearity well by sharing weight, always has a unique solution, but never yields exact zeros. Lasso (L1) yields sparse models via soft-thresholding, useful when few features matter, but picks arbitrarily among correlated features and selects at most n features when p > n. Elastic Net mixes both: sparse and stable with correlated groups. Choose by the goal (interpretability/sparsity vs pure prediction) and by cross-validated error; standardise features first.

Q8. Why is polynomial regression still "linear" regression?
Linearity refers to the parameters. y = β₀ + β₁x + β₂x² is a linear combination of the columns 1, x, x² with coefficients β, so OLS, the normal equations and the usual inference apply. The fitted curve is non-linear in x. Caveats: powers are collinear (center or use orthogonal polynomials), high degree overfits and extrapolates badly.

Q9. How do you avoid data leakage in a regression workflow?
Split first (or use CV), then learn every data-dependent transformation (imputation, scaling, encoding, feature selection, target encoding) inside a Pipeline so it is fit only on training folds. Use TimeSeriesSplit for temporal data and GroupKFold when the same entity appears in many rows. Exclude features that would not be available at prediction time. Use the test set once.

Q10. A linear regression model is in production. What do you monitor and when do you retrain?
Inputs: schema violations, missing rates, out-of-range rates, feature distribution drift (for example PSI or KS tests against the training distribution). Outputs: prediction distribution. Performance once labels arrive: rolling RMSE/MAE, residual bias overall and by segment, prediction-interval coverage. Operations: latency and errors. Retrain on a schedule or when drift/error thresholds are crossed, validate the challenger against the current model on recent data, deploy gradually (shadow or canary), and keep rollback. Version the entire pipeline, data snapshot and library versions.

4. Mini project brief — "Explainable price model with honest uncertainty"

Dataset: California housing (sklearn.datasets.fetch_california_housing). It is downloaded once from the internet on first use (then cached under ~/scikit_learn_data); everything else runs locally at no cost. In a check run it returned 20,640 rows, 8 numeric features (MedInc, HouseAge, AveRooms, AveBedrms, Population, AveOccup, Latitude, Longitude) and target MedHouseVal (median house value in units of $100,000), with a maximum of about 5.0, which suggests the target is capped. If network access is not available, use the diabetes dataset from Lab 3.

Deliverables 1. A reproducible script or notebook with fixed random seeds and a pinned requirements.txt. 2. An EDA section: distributions, correlation matrix, VIF table; note the capped target and decide (with justification) how to treat capped rows. 3. A baseline Pipeline(StandardScaler, LinearRegression) with 5-fold CV (report mean ± SD of RMSE and R²) and a single held-out test evaluation. 4. At least two feature-engineering iterations (for example log of skewed features, AveBedrms/AveRooms ratio, polynomial or interaction terms via PolynomialFeatures, spatial bins of latitude/longitude) with CV evidence for each. 5. RidgeCV, LassoCV and ElasticNetCV comparison inside pipelines; a coefficient path or coefficient table on the standardised scale. 6. statsmodels OLS of the final feature set: summary, Breusch–Pagan, Q-Q plot, Cook's distance, and a comparison of classical vs HC3 standard errors. 7. Prediction intervals for 5 example rows and an empirical coverage check of 95% prediction intervals on the test set. 8. A one-page model card: intended use, data, metrics, limitations (capping, spatial correlation, extrapolation), monitoring plan and retraining trigger.

Assessment rubric (suggested) - Correct, leakage-free evaluation (30%) - Diagnostics and correct statistical interpretation (25%) - Regularisation and feature engineering backed by CV evidence (20%) - Production readiness: pipeline serialisation, input validation, monitoring plan (15%) - Clarity of the model card (10%)

5. Flashcards

#	Question	Answer
1	Normal equations?	XᵀXβ̂ = Xᵀy
2	Trace of the hat matrix?	p + 1, the number of fitted parameters
3	Two properties of OLS residuals (with intercept)?	Sum to zero; orthogonal to every column of X (Xᵀe = 0)
4	Relationship between cond(X) and cond(XᵀX)?	cond(XᵀX) = cond(X)²
5	What does BLUE stand for?	Best Linear Unbiased Estimator
6	Is normality a Gauss–Markov assumption?	No; it is needed only for exact small-sample t/F distributions
7	Effect of heteroscedasticity on OLS coefficients?	None on bias; it invalidates classical standard errors and removes efficiency
8	VIF formula and SE inflation?	VIF_j = 1/(1 − R²_j); SE inflated by √VIF_j
9	Durbin–Watson value for no autocorrelation?	About 2 (range 0–4; DW ≈ 2(1 − ρ̂))
10	Cook's distance rules of thumb?	> 1 large; > 4/n worth inspecting
11	Adjusted R² formula?	1 − (1 − R²)(n − 1)/(n − p − 1)
12	Ridge closed form?	(XᵀX + λI)⁻¹Xᵀy (centered data, intercept unpenalised)
13	Why does lasso give exact zeros?	L1 penalty is non-differentiable at 0 → soft-thresholding; L1 ball has corners on the axes
14	Fix for the dummy variable trap?	Use k − 1 dummies with an intercept (or drop the intercept)
15	Confidence vs prediction interval?	CI is for the mean response; PI is for a new observation and adds σ², so it is wider

6. Further reading (official documentation and standard texts)

All links were checked with curl and returned HTTP 200 on 2026-10-02.

scikit-learn - Linear models user guide (OLS, Ridge, Lasso, Elastic Net): https://scikit-learn.org/stable/modules/linear_model.html - LinearRegression: https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html - RidgeCV: https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.RidgeCV.html - LassoCV: https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LassoCV.html - ElasticNetCV: https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNetCV.html - Pipeline: https://scikit-learn.org/stable/modules/generated/sklearn.pipeline.Pipeline.html - Cross-validation: https://scikit-learn.org/stable/modules/cross_validation.html - Common pitfalls (including data leakage): https://scikit-learn.org/stable/common_pitfalls.html - Toy datasets (diabetes): https://scikit-learn.org/stable/datasets/toy_dataset.html

statsmodels - OLS: https://www.statsmodels.org/stable/generated/statsmodels.regression.linear_model.OLS.html - Linear regression overview: https://www.statsmodels.org/stable/regression.html - variance_inflation_factor: https://www.statsmodels.org/stable/generated/statsmodels.stats.outliers_influence.variance_inflation_factor.html - het_breuschpagan: https://www.statsmodels.org/stable/generated/statsmodels.stats.diagnostic.het_breuschpagan.html - durbin_watson: https://www.statsmodels.org/stable/generated/statsmodels.stats.stattools.durbin_watson.html - Regression diagnostics example: https://www.statsmodels.org/stable/examples/notebooks/generated/regression_diagnostics.html

NumPy - numpy.linalg.lstsq: https://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html - numpy.linalg.pinv: https://numpy.org/doc/stable/reference/generated/numpy.linalg.pinv.html - numpy.linalg.qr: https://numpy.org/doc/stable/reference/generated/numpy.linalg.qr.html

Books (free online editions from the authors) - James, Witten, Hastie, Tibshirani (and Taylor), An Introduction to Statistical Learning: https://www.statlearning.com/ - Hastie, Tibshirani, Friedman, The Elements of Statistical Learning: https://hastie.su.domains/ElemStatLearn/