🎯 What You'll Learn
- Understand regression as the task of predicting a continuous numerical output
- Derive the Ordinary Least Squares (OLS) solution mathematically and implement it with NumPy
- Understand the difference between MSE, MAE, and RMSE, and when to choose each
- Interpret R² and Adjusted R² to evaluate how much variance a model explains
- Train, evaluate, and interpret a linear regression model using scikit-learn
pip install scikit-learn numpy matplotlib pandas. This lesson uses mathematical notation: ŷ (y-hat) means predicted value; y means actual value; β (beta) means model parameter. These conventions are standard across all ML literature.
1 What Is Regression?
Regression is a type of supervised learning where the target variable is continuous and numerical — not a category. The model's job is to learn a mapping from input features to a real-valued output.
Examples of regression problems:
- House price prediction: features (sqft, bedrooms, location) → output (sale price in $)
- Salary estimation: features (years experience, education, role) → output (annual salary in $)
- Weather forecasting: features (pressure, humidity, temperature yesterday) → output (temperature tomorrow in °C)
- Stock returns: features (technical indicators, volume, sentiment) → output (next-day return %)
This is distinct from classification, where the output is a discrete category. If you're predicting "will this person default (yes/no)?", that's classification. If you're predicting "what is the probability of default?", that's still classification (0–1 bounded output). If you're predicting "what will the default rate be across 10,000 loans?", that's regression.
Regression: continuous output (price, temperature, score) — measured with the error metrics you'll learn in this lesson (RMSE, MAE, R²). Classification: discrete output (class label, probability) — measured with its own family of metrics you'll meet in Lessons 19–20. The boundary blurs when you predict a probability (0.73), but the task — and therefore the algorithm and loss function — depends on whether you are predicting a quantity or a category.
2 Simple Linear Regression
The simplest regression model assumes a linear relationship between one input feature x and the output y:
ŷ = β₀ + β₁x
where:
- β₀ (intercept) — the predicted value of y when x = 0; shifts the line up or down
- β₁ (slope) — how much y changes for a one-unit increase in x; the "effect size" of the feature
- ŷ — the model's prediction (read: "y-hat")
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
# Generate synthetic data: salary vs years_experience
np.random.seed(42)
n = 100
years_experience = np.random.uniform(0, 15, n)
salary = 30000 + 4500 * years_experience + np.random.normal(0, 8000, n)
X = years_experience.reshape(-1, 1) # must be 2D for sklearn
y = salary
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Fit simple linear regression
model = LinearRegression()
model.fit(X_train, y_train)
print(f"Intercept (β₀): {model.intercept_:.2f}") # ≈ 30,000
print(f"Slope (β₁): {model.coef_[0]:.2f}") # ≈ 4,500
# Predict
y_pred = model.predict(X_test)
print(f"Prediction for 5 years exp: ${model.predict([[5]])[0]:,.0f}")
# ≈ $52,500 (30,000 + 4,500 × 5)The slope β₁ ≈ 4,500 means: each additional year of experience is associated with a $4,500 increase in predicted salary, holding all else constant. This direct interpretability is one of linear regression's greatest strengths.
3 Ordinary Least Squares (OLS): Finding the Best Line
What makes a line the "best fit"? For each training point, the model makes a prediction ŷᵢ. The difference between the actual value yᵢ and the prediction is called the residual:
residualᵢ = yᵢ − ŷᵢ
OLS minimises the sum of squared residuals (also called the Residual Sum of Squares, RSS):
RSS = Σᵢ (yᵢ − ŷᵢ)² = Σᵢ (yᵢ − β₀ − β₁xᵢ)²
By taking the partial derivatives of RSS with respect to β₀ and β₁ and setting them to zero, we get the closed-form (analytical) solution. In matrix notation for multiple features:
β = (XᵀX)⁻¹ Xᵀy
where X is the feature matrix (with a bias column of ones prepended), y is the target vector, and β is the vector of coefficients.
Try It Yourself: Fit a Line by Hand
Before trusting the formula, build intuition for what it's actually doing. Below is a small salary dataset (years of experience → salary, the same setup as Section 2). Drag the slope and intercept sliders to draw your own candidate line through the points, and watch the MSE update live. Try to get the MSE as low as you can manually — then compare to the OLS-optimal line marked on the chart.
Your line: ŷ = 30,000 + 4,000x — drag the sliders and watch the MSE change.
The dashed green line marks the OLS-optimal fit (computed once via the normal equation) — the unique line that minimises MSE over this dataset. No matter how you drag the sliders, you cannot beat its MSE; that's exactly what "OLS minimises RSS" means in practice.
Why Squared Residuals (not absolute)?
- Differentiable everywhere: the squared function is smooth — we can take its derivative at every point. The absolute value function has a sharp corner at 0 and is not differentiable there. This matters for gradient-based optimization.
- Penalises large errors more: squaring amplifies large errors. If a residual is 10 instead of 1, squared error is 100× larger (100 vs 1), but absolute error is only 10× larger. This makes OLS focus on reducing large mistakes.
- Unique closed-form solution: minimizing the sum of squares has one exact mathematical answer (assuming the columns of X are linearly independent). Minimizing sum of absolute values does not.
import numpy as np
def ols_from_scratch(X, y):
"""
Compute OLS solution: beta = (X^T X)^{-1} X^T y
X: feature matrix, shape (n_samples, n_features)
y: target vector, shape (n_samples,)
"""
# Add bias column (column of 1s) as first column
n = X.shape[0]
X_b = np.column_stack([np.ones(n), X]) # shape: (n, n_features + 1)
# Normal equation: (X^T X)^{-1} X^T y
beta = np.linalg.inv(X_b.T @ X_b) @ X_b.T @ y
return beta # beta[0] = intercept, beta[1:] = coefficients
# Test on our salary data
X_arr = years_experience.reshape(-1, 1)
beta = ols_from_scratch(X_arr, salary)
print(f"OLS Intercept: {beta[0]:.2f}") # ≈ 30,000
print(f"OLS Slope: {beta[1]:.2f}") # ≈ 4,500
# Verify: matches sklearn exactly
from sklearn.linear_model import LinearRegression
lr = LinearRegression()
lr.fit(X_arr, salary)
print(f"sklearn: intercept={lr.intercept_:.2f}, coef={lr.coef_[0]:.2f}")The normal equation requires computing (XᵀX)⁻¹, which fails if the feature matrix is singular (non-invertible). This happens when features are perfectly multicollinear — when one feature is a linear combination of others. For example, if you include both "income" and "income_in_thousands" in the same model. Fix: remove redundant features, or use sklearn, which handles this gracefully (Lesson 21 adds a third option, regularization). For very large datasets, OLS is also too slow (O(n·p²) complexity) — the fix for that is the subject of the very next lesson: gradient descent.
4 Cost Functions: MSE, MAE, and RMSE
A cost function measures how wrong your model's predictions are. You compute it on a set of predictions and compare it to actual values. The choice of cost function affects what you optimize for during training.
Mean Squared Error (MSE)
MSE = (1/n) Σᵢ (yᵢ − ŷᵢ)²
MSE is the average squared error. It is sensitive to outliers because errors are squared — a single prediction that's off by 100 contributes 10,000 to MSE. The units are the squared units of your target variable (e.g., $² for salary), which is hard to interpret directly.
Root Mean Squared Error (RMSE)
RMSE = √MSE = √[ (1/n) Σᵢ (yᵢ − ŷᵢ)² ]
RMSE restores the original units by taking the square root of MSE. An RMSE of $8,000 means, on average, predictions are off by about $8,000. RMSE is the most commonly reported regression metric because it's in the same unit as the target.
Mean Absolute Error (MAE)
MAE = (1/n) Σᵢ |yᵢ − ŷᵢ|
MAE is the average of the absolute differences between predictions and actual values. It treats all errors linearly — an error of 100 contributes exactly 100×, not 10,000×. This makes MAE more robust to outliers. An MAE of $5,000 means the average prediction is off by $5,000.
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
import numpy as np
# Actual vs predicted salaries
y_true = np.array([45000, 60000, 75000, 50000, 90000])
y_pred = np.array([42000, 65000, 71000, 53000, 88000])
mse = mean_squared_error(y_true, y_pred)
rmse = mean_squared_error(y_true, y_pred, squared=False) # or np.sqrt(mse)
mae = mean_absolute_error(y_true, y_pred)
print(f"MSE: {mse:,.0f}") # 19,200,000 (in $² — hard to interpret)
print(f"RMSE: {rmse:,.0f}") # ~4,382 (in $ — "on average off by $4,382")
print(f"MAE: {mae:,.0f}") # 3,400 (in $ — "average absolute error $3,400")
# Manual verification
errors = y_true - y_pred # [-3000, 5000, -4000, 3000, -2000]
print(np.mean(errors ** 2)) # MSE
print(np.sqrt(np.mean(errors ** 2))) # RMSE
print(np.mean(np.abs(errors))) # MAE| Metric | Outlier sensitivity | Units | Use when |
|---|---|---|---|
| MSE | High (penalises large errors heavily) | Squared (y²) | Training loss during optimization (next lesson) |
| RMSE | High (same as MSE) | Same as y ⭐ | Reporting model performance (interpretable) |
| MAE | Low (robust to outliers) | Same as y | Data with outliers; business-facing reports |
5 Model Evaluation: R² and Adjusted R²
MSE and RMSE are absolute — they tell you the magnitude of error, but they don't tell you if the model is actually good relative to simply predicting the mean. R² (R-squared, or the coefficient of determination) answers this question.
R² = 1 − SS_res / SS_tot where SS_res = Σ(y − ŷ)² SS_tot = Σ(y − ȳ)²
Interpretation:
- R² = 1.0: perfect predictions — all variance is explained
- R² = 0.0: the model does no better than simply predicting the mean ȳ for every point
- R² = 0.75: the model explains 75% of the variance in y; 25% remains unexplained
- R² < 0: the model is worse than predicting the mean (this can happen in extreme cases or on the test set)
from sklearn.metrics import r2_score
import numpy as np
y_true = np.array([45000, 60000, 75000, 50000, 90000, 55000])
y_pred = np.array([42000, 63000, 73000, 51000, 87000, 57000])
r2 = r2_score(y_true, y_pred)
print(f"R²: {r2:.4f}") # e.g., 0.9821
# Manual calculation
y_mean = y_true.mean()
ss_res = np.sum((y_true - y_pred) ** 2) # residual sum of squares
ss_tot = np.sum((y_true - y_mean) ** 2) # total sum of squares
r2_manual = 1 - ss_res / ss_tot
print(f"R² (manual): {r2_manual:.4f}")Adjusted R²: Penalizing Unnecessary Features
A problem with R²: it never decreases when you add features, even if those features are useless noise. Adjusted R² corrects for this by penalizing the addition of features that don't meaningfully improve fit:
Adjusted R² = 1 − (1 − R²) × (n − 1) / (n − p − 1)
where n = number of samples, p = number of features. If a new feature improves the model meaningfully, Adjusted R² increases. If it's noise, Adjusted R² decreases or stays flat. Use Adjusted R² when comparing models with different numbers of features.
def adjusted_r2(r2, n, p):
"""
r2: regular R-squared
n: number of data points
p: number of features (not counting intercept)
"""
return 1 - (1 - r2) * (n - 1) / (n - p - 1)
n, p = 100, 5
r2 = 0.82
adj_r2 = adjusted_r2(r2, n, p)
print(f"R²: {r2:.4f}")
print(f"Adjusted R²: {adj_r2:.4f}") # slightly lower because we have 5 features
# If we add a useless noise feature (p=6), R² stays the same but adj_r2 drops
r2_noise = 0.8201 # barely improves
adj_r2_noise = adjusted_r2(r2_noise, n, p=6)
print(f"Adj R² with noise feature: {adj_r2_noise:.4f}") # lower!6 Implementing with scikit-learn
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
# Generate multivariate salary dataset
np.random.seed(42)
n = 500
data = {
'years_experience': np.random.uniform(0, 20, n),
'education_years': np.random.uniform(12, 22, n),
'num_skills': np.random.randint(1, 15, n),
}
df = pd.DataFrame(data)
# True relationship + noise
salary = (25000
+ 4500 * df['years_experience']
+ 2000 * df['education_years']
+ 800 * df['num_skills']
+ np.random.normal(0, 5000, n))
df['salary'] = salary
X = df[['years_experience', 'education_years', 'num_skills']]
y = df['salary']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Pipeline: scale + linear regression
pipeline = Pipeline([
('scaler', StandardScaler()),
('model', LinearRegression()),
])
pipeline.fit(X_train, y_train)
y_pred = pipeline.predict(X_test)
# Evaluation
rmse = mean_squared_error(y_test, y_pred, squared=False)
mae = mean_absolute_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
print(f"RMSE: ${rmse:,.0f}") # ≈ $5,100
print(f"MAE: ${mae:,.0f}") # ≈ $4,100
print(f"R²: {r2:.4f}") # ≈ 0.91
# Access coefficients (note: scaling changes coef magnitude in scaled space)
lr = pipeline.named_steps['model']
scaler = pipeline.named_steps['scaler']
print("\nCoefficients (in original feature units):")
for name, coef, scale in zip(X.columns, lr.coef_, scaler.scale_):
coef_original = coef / scale # unscale
print(f" {name}: {coef_original:.1f}")When you apply StandardScaler before LinearRegression, the model.coef_ values are in standardized units (i.e., "per standard deviation change in that feature"). To recover the original-unit coefficients (e.g., "per year of experience"), divide each coefficient by the corresponding scaler.scale_ value. The direct sklearn coefficients are useful for comparing relative importance of features — the feature with the largest |coef| after scaling contributes most to predictions.
7 Assumptions of Linear Regression
OLS gives optimal estimates only when five assumptions hold. Violations don't make the algorithm "wrong", but they mean the coefficients may be biased or the confidence intervals unreliable.
- Linearity: the relationship between X and y is linear. Check with residual plot — if you see a U-shaped pattern, linearity is violated.
- Independence of errors: residuals for different observations are not correlated. Especially important for time-series data.
- Homoscedasticity: the variance of residuals is constant across all fitted values. If residuals fan out as ŷ increases, the assumption is violated.
- Normality of residuals: residuals are approximately normally distributed. Required for valid p-values and confidence intervals, not for predictions.
- No perfect multicollinearity: no feature is a perfect linear combination of other features. Imperfect multicollinearity is fine in practice.
Here's the residuals-vs-fitted diagnostic applied to the OLS line from the slider exercise above — fitted on the same 25-point salary dataset:
Residuals vs fitted values for the OLS-optimal line. Points scatter randomly around the zero line with no obvious funnel or curve — consistent with linearity and homoscedasticity holding reasonably well.
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression
# After fitting: check residuals
np.random.seed(42)
X_data = np.random.randn(200, 1)
y_data = 3 * X_data.ravel() + np.random.randn(200) * 2 # linear + noise
model = LinearRegression().fit(X_data, y_data)
y_hat = model.predict(X_data)
residuals = y_data - y_hat
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
# Residuals vs fitted values (check linearity + homoscedasticity)
axes[0].scatter(y_hat, residuals, alpha=0.5, color='steelblue')
axes[0].axhline(0, color='red', linestyle='--')
axes[0].set_xlabel('Fitted values')
axes[0].set_ylabel('Residuals')
axes[0].set_title('Residuals vs Fitted')
# Histogram of residuals (check normality)
axes[1].hist(residuals, bins=20, color='steelblue', edgecolor='white')
axes[1].set_xlabel('Residual')
axes[1].set_title('Distribution of Residuals')
plt.tight_layout()
plt.savefig('residual_plots.png', dpi=100)
print("Saved to residual_plots.png")
# A random scatter around 0 in plot 1 = assumptions satisfied
# A bell curve in plot 2 = normality satisfied8 From Scratch with NumPy: The Normal Equation
Implementing the OLS solution from scratch solidifies your understanding and makes the connection between the math and the code explicit:
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.datasets import make_regression
# Generate regression data
X, y = make_regression(n_samples=200, n_features=3, noise=15, random_state=42)
class LinearRegressionOLS:
"""Linear Regression via the Normal Equation (closed-form OLS)."""
def fit(self, X, y):
# Add intercept column (column of 1s at position 0)
n = X.shape[0]
X_b = np.column_stack([np.ones(n), X]) # shape: (n, p+1)
# Normal equation: beta = (X^T X)^{-1} X^T y
# Use lstsq instead of direct inv() for numerical stability
self.beta_, residuals, rank, sv = np.linalg.lstsq(X_b, y, rcond=None)
self.intercept_ = self.beta_[0]
self.coef_ = self.beta_[1:]
return self
def predict(self, X):
return X @ self.coef_ + self.intercept_
def score(self, X, y):
y_pred = self.predict(X)
ss_res = np.sum((y - y_pred) ** 2)
ss_tot = np.sum((y - y.mean()) ** 2)
return 1 - ss_res / ss_tot
# Compare to sklearn
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
ols_custom = LinearRegressionOLS().fit(X_train, y_train)
ols_sklearn = LinearRegression().fit(X_train, y_train)
print("Custom OLS R²:", ols_custom.score(X_test, y_test))
print("sklearn LR R²:", ols_sklearn.score(X_test, y_test))
# Both should match to many decimal places
print("\nCoefficients match?")
print("Custom: ", ols_custom.coef_.round(4))
print("sklearn:", ols_sklearn.coef_.round(4))Real-World Spotlight: Employee Salary Prediction
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler, OrdinalEncoder, OneHotEncoder
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_absolute_error, r2_score
# Build dataset
np.random.seed(42)
n = 300
education_map = {'BS': 0, 'MS': 1, 'PhD': 2}
edu_choices = np.random.choice(['BS', 'MS', 'PhD'], n, p=[0.5, 0.35, 0.15])
dept_choices = np.random.choice(['Engineering', 'Sales', 'Marketing'], n, p=[0.6, 0.25, 0.15])
df = pd.DataFrame({
'years_exp': np.random.uniform(0, 20, n),
'education': edu_choices,
'department': dept_choices,
})
# True salary formula
edu_bonus = np.array([education_map[e] for e in edu_choices]) * 8000
dept_bonus = np.where(df['department'] == 'Engineering', 15000,
np.where(df['department'] == 'Sales', 5000, 0))
df['salary'] = (35000
+ 4200 * df['years_exp']
+ edu_bonus
+ dept_bonus
+ np.random.normal(0, 6000, n))
X = df[['years_exp', 'education', 'department']]
y = df['salary']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Build preprocessing + model pipeline
preprocessor = ColumnTransformer([
('num', StandardScaler(), ['years_exp']),
('ord', OrdinalEncoder(categories=[['BS', 'MS', 'PhD']]), ['education']),
('cat', OneHotEncoder(drop='first', sparse_output=False), ['department']),
])
pipeline = Pipeline([
('prep', preprocessor),
('model', LinearRegression()),
])
pipeline.fit(X_train, y_train)
y_pred = pipeline.predict(X_test)
mae = mean_absolute_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
print(f"MAE: ${mae:,.0f}") # ≈ $4,800
print(f"R²: {r2:.4f}") # ≈ 0.90
# Interpret: how much does a PhD add over a BS?
# Coef for 'education' feature × OrdinalEncoder spacing = per-level value
lr = pipeline.named_steps['model']
feature_names = ['years_exp_scaled', 'education_encoded', 'dept_Sales', 'dept_Marketing']
print("\nCoefficients:")
for name, coef in zip(feature_names, lr.coef_):
print(f" {name}: {coef:.1f}")Quick Check
✍️ Practice Exercises
- Load the California Housing dataset (
from sklearn.datasets import fetch_california_housing). Train a LinearRegression model using a Pipeline with StandardScaler. Report RMSE, MAE, and R² on the test set. Which feature has the largest coefficient? - Implement the OLS normal equation from scratch using NumPy (without using LinearRegression). Apply it to a dataset with 2 features. Verify that your intercept and coefficients match sklearn's LinearRegression within 4 decimal places.
- Generate 50 data points with a quadratic relationship (y = 2x² + noise). Fit a LinearRegression model. Plot the residuals. What pattern do you see, and what does it tell you about the linear assumption?
- Compare MSE vs MAE on a dataset where you introduce one outlier (multiply one y_true value by 10). Measure how each metric changes. Which is more affected? Which should you use in that scenario?
▶ Show California Housing solution
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
import numpy as np
data = fetch_california_housing()
X, y = data.data, data.target
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
pipeline = Pipeline([('scaler', StandardScaler()), ('model', LinearRegression())])
pipeline.fit(X_train, y_train)
y_pred = pipeline.predict(X_test)
rmse = mean_squared_error(y_test, y_pred, squared=False)
mae = mean_absolute_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
print(f"RMSE: {rmse:.4f}, MAE: {mae:.4f}, R²: {r2:.4f}")
# Feature importances
lr = pipeline.named_steps['model']
for name, coef in sorted(zip(data.feature_names, lr.coef_), key=lambda x: abs(x[1]), reverse=True):
print(f" {name}: {coef:.4f}")📚 Primary Source for This Lesson
scikit-learn User Guide: Linear Models
Covers LinearRegression, Ridge, Lasso, and other linear models with mathematical formulations and usage examples. Also see Chapter 3 of An Introduction to Statistical Learning (freely available at statlearning.com) — the gold standard reference for regression and its assumptions.