Linear Regression: From Least Squares to Production-Ready Practice

Linear Regression: From Least Squares to Production-Ready Practice

Tags: machinelearning, datascience, python, tutorial

Linear regression is the first algorithm most people learn, and the one most people never study deeply. It is also the model you will still find in production after fancier algorithms fail, because it is fast, stable, and explainable.

This article is not a "call .fit() and read the score" tutorial. We will cover the math, the statistical assumptions, the diagnostics, regularization, evaluation, production concerns, and the interview questions that separate beginners from engineers.

Why Linear Regression Deserves a Second Look

Linear regression is the foundation for understanding almost every other supervised model:

• Logistic regression is linear regression with a sigmoid on top.
• Ridge and Lasso are linear regression with constrained weights.
• Neural networks are stacked linear transformations with nonlinear activations.
• Tree models are judged against the same baseline: "can I beat a linear model?"

More importantly, linear regression is still the right answer in many business problems. When you need to explain a prediction to a regulator, a client, or a finance team, a clean linear model with interpretable coefficients beats a black box.

The Math: Least Squares and the Normal Equation

Given features X and target y, a linear model assumes:

y = X * beta + epsilon

Enter fullscreen mode Exit fullscreen mode

The goal is to minimize the residual sum of squares:

L(beta) = ||y - X*beta||^2

Enter fullscreen mode Exit fullscreen mode

Taking the derivative with respect to beta and setting it to zero gives the normal equation:

beta = (X^T * X)^(-1) * X^T * y

Enter fullscreen mode Exit fullscreen mode

In practice, use the pseudoinverse (pinv) instead of the inverse, because X^T X may be singular or numerically unstable when features are collinear.

import numpy as np

def normal_equation(X, y):
    Xb = np.c_[np.ones(X.shape[0]), X]  # add intercept
    beta = np.linalg.pinv(Xb.T @ Xb) @ Xb.T @ y
    return beta

Enter fullscreen mode Exit fullscreen mode

Three Equivalent Views of Least Squares

1. Geometric view

The prediction X * beta is a projection of y onto the column space of X. Least squares finds the closest point in that subspace. This is why residuals are orthogonal to the fitted values.

2. Maximum likelihood view

Assume:

epsilon ~ N(0, sigma^2)

Enter fullscreen mode Exit fullscreen mode

Then maximizing the log-likelihood is equivalent to minimizing the sum of squared errors. This connection explains why normality of residuals matters for confidence intervals, even though OLS coefficients are still consistent under milder assumptions.

3. Optimization view

For large datasets, computing (X^T X)^(-1) becomes expensive. Gradient descent solves the same objective iteratively:

def gradient_descent(X, y, lr=0.01, epochs=500):
    Xb = np.c_[np.ones(X.shape[0]), X]
    n, d = Xb.shape
    beta = np.zeros(d)
    for _ in range(epochs):
        grad = (2 / n) * Xb.T @ (Xb @ beta - y)
        beta -= lr * grad
    return beta

Enter fullscreen mode Exit fullscreen mode

Normal equation: exact, best for small and medium datasets. Gradient descent: scalable, but needs feature scaling and learning-rate tuning.

Assumptions You Must Check

A fitted model with a high R-squared is not automatically trustworthy. These assumptions are what make the coefficients interpretable:

Assumption	Meaning	How to check
Linearity	Linear relationship between features and target	residual vs fitted plot
Independence	Errors are not correlated	Durbin-Watson statistic
Homoscedasticity	Constant error variance	Breusch-Pagan test, residual plot
Normality	Errors are normally distributed	Q-Q plot, Jarque-Bera test
No multicollinearity	Features are not highly correlated	VIF, correlation matrix

Complete diagnostic walkthrough

import statsmodels.api as sm
from statsmodels.stats.stattools import durbin_watson
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.outliers_influence import variance_inflation_factor
from scipy import stats
import pandas as pd

X_const = sm.add_constant(X_train)
model = sm.OLS(y_train, X_const).fit()
resid = model.resid

print("Durbin-Watson:", durbin_watson(resid))
# Values near 2 mean no autocorrelation; near 0 or 4 is a warning.

lm, lm_pvalue, fvalue, f_pvalue = het_breuschpagan(resid, X_const)
print("Breusch-Pagan p-value:", f_pvalue)
# p > 0.05: no strong evidence of heteroscedasticity.

print("Jarque-Bera p-value:", stats.jarque_bera(resid).pvalue)
# p > 0.05: no strong evidence against normality.

vif = pd.DataFrame({
    "feature": X_train.columns,
    "VIF": [variance_inflation_factor(X_train.values, i)
            for i in range(X_train.shape[1])],
})
print(vif)
# VIF > 10 is commonly treated as a multicollinearity warning.

Enter fullscreen mode Exit fullscreen mode

Regularization: The Bias-Variance Tradeoff

Plain least squares minimizes training error and can overfit when features are many or collinear. Regularization adds a penalty:

Ridge:       L = MSE + alpha * sum(beta_i^2)
Lasso:       L = MSE + alpha * sum(|beta_i|)
ElasticNet:  L = MSE + alpha * (r * L1 + (1 - r) * L2)

Enter fullscreen mode Exit fullscreen mode

Model	Penalty	Effect	Best when
Ridge	L2	Shrinks weights, keeps all features	Many correlated features
Lasso	L1	Shrinks some weights to exactly zero	Feature selection needed
ElasticNet	L1 + L2	Sparse but stable with correlated groups	Mix of both needs

Critical detail: always standardize features before regularization. Otherwise the penalty unfairly shrinks large-magnitude features.

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet

models = {
    "OLS": LinearRegression(),
    "Ridge(alpha=1)": Ridge(alpha=1.0),
    "Lasso(alpha=0.01)": Lasso(alpha=0.01),
    "ElasticNet(alpha=0.01)": ElasticNet(alpha=0.01, l1_ratio=0.5),
}

for name, model in models.items():
    pipe = make_pipeline(StandardScaler(), model)
    pipe.fit(X_train, y_train)
    y_pred = pipe.predict(X_test)
    print(f"{name:20s} RMSE={mean_squared_error(y_test, y_pred, squared=False):.4f} "
          f"R2={r2_score(y_test, y_pred):.4f}")

Enter fullscreen mode Exit fullscreen mode

Evaluation Metrics

Do not report only R-squared. Each metric tells a different story:

Metric	Formula meaning	Use when
MAE	Mean absolute error	Errors should not be squared, outliers matter less
RMSE	Root mean squared error	Large errors are especially bad
MAPE	Mean absolute percentage error	Business wants a percentage interpretation
R²	Proportion of variance explained	Comparing model quality
Adjusted R²	R² penalized by feature count	Comparing models with different feature sets

Always compare against a mean baseline. If your model only beats "predict the average" by a little, the features are not adding much.

baseline_pred = np.full_like(y_test, y_train.mean())
print("Baseline RMSE:", mean_squared_error(y_test, baseline_pred, squared=False))

Enter fullscreen mode Exit fullscreen mode

Complete End-to-End Example

import numpy as np
import pandas as pd
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

data = fetch_california_housing(as_frame=True)
X = data.frame.drop(columns=["MedHouseVal"])
y = data.frame["MedHouseVal"]

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# 1. Handcrafted normal equation
beta = normal_equation(X_train.values, y_train.values)
X_test_b = np.c_[np.ones(X_test.shape[0]), X_test.values]
y_pred_manual = X_test_b @ beta

# 2. sklearn baseline
from sklearn.linear_model import LinearRegression
lr = LinearRegression().fit(X_train, y_train)
y_pred_sklearn = lr.predict(X_test)

for name, y_pred in [("manual", y_pred_manual), ("sklearn", y_pred_sklearn)]:
    print(f"{name}: RMSE={mean_squared_error(y_test, y_pred, squared=False):.4f} "
          f"MAE={mean_absolute_error(y_test, y_pred):.4f} "
          f"R2={r2_score(y_test, y_pred):.4f}")

Enter fullscreen mode Exit fullscreen mode

The manual solution and the sklearn solution should produce nearly identical numbers. If they do not, your gradient descent is not converged or your features are not scaled.

Production Checklist

1. Feature scaling: required for gradient descent and regularized models.
2. Cross-validation: never tune alpha on the test set.
3. Data leakage: split before scaling, encoding, or imputation.
4. Monitoring: track feature distributions and prediction drift after deployment.
5. Explainability: record coefficients, feature importance, and SHAP values for audits.
6. Simplicity first: start with linear regression before reaching for XGBoost.

Common Mistakes

1. Correlation is not causation

A high coefficient does not mean changing the feature will change the target. Omitted variables and reverse causality are always possible.

2. Extrapolation

Linear models are dangerous outside the training range. A model trained on houses up to 300 sqm will not tell you what a 3000 sqm house costs.

3. High R-squared does not mean good prediction

R-squared measures in-sample fit. A model can have high R-squared on training data and terrible RMSE on unseen data.

4. Ignoring residuals

If residuals have a pattern, your model is missing structure. Do not "fix" a nonlinear pattern by adding more linear features blindly.

5. Using raw features with regularization

Without standardization, Ridge and Lasso apply unequal penalties to features measured in different units.

Interview Questions to Be Ready For

1. What is the difference between the normal equation and gradient descent?
2. Why is OLS equivalent to maximum likelihood under a normal error assumption?
3. What happens if X^T X is singular, and how does Ridge help?
4. How do you detect and fix multicollinearity?
5. Why use adjusted R-squared instead of R-squared?
6. What does a residual vs fitted plot tell you?
7. When would you prefer Lasso over Ridge?
8. Why must you standardize before L1 or L2 regularization?
9. What are the consequences of heteroscedasticity?
10. How do you detect data leakage in a regression pipeline?

Conclusion

Linear regression looks simple, but it contains the core ideas of machine learning: objective functions, optimization, statistical assumptions, regularization, evaluation, and deployment. If you deeply understand this one model, every algorithm after it becomes easier to learn.

Start with the math, validate the assumptions, and keep the production concerns in mind. The model that is "too simple" is often the one that survives in production the longest.

---

Found this useful? Follow me for more practical AI and data engineering posts.
没 ；。l