Ordinary Least Squares (OLS) Regression: The Complete Guide

1. Introduction

Regression analysis is the statistical practice of asking a single, disciplined question: how does one number change when other numbers change? When a marketing analyst wonders how much revenue climbs for every extra dollar spent on ads, when a hospital models how blood pressure moves with age, or when an insurer estimates how medical charges scale with body mass index — each of them is running, formally or informally, a regression. Regression is the vocabulary we use to describe relationships between a target variable and one or more explanatory variables.

Ordinary Least Squares — usually abbreviated OLS — is the oldest, most studied, and still most widely used technique for fitting a linear regression model. The name tells you exactly what it does: it chooses the straight line (or hyperplane, in higher dimensions) that makes the sum of the squared vertical distances between the observed data points and the line as small as possible. That is the entire idea. Everything else — assumptions, diagnostics, p-values, confidence intervals — is machinery built on top of that single geometric commitment.

OLS matters far beyond its own accuracy scoreboard. It is the intellectual scaffolding for logistic regression, generalized linear models, mixed-effects models, Ridge, Lasso, Elastic Net, and much of classical econometrics. If you understand OLS deeply, most of the rest of statistical learning becomes an exercise in bending the same ideas to fit a new shape.

Prediction vs. Statistical Inference

It is worth being explicit about a subtle distinction. Machine learning often cares only about prediction — given a new input, produce the most accurate output. Classical statistics cares about inference — understanding which inputs actually matter, in what direction, by how much, and with what uncertainty. OLS is a rare technique that speaks both languages fluently. A single fitted model gives you both a prediction machine and a rich set of hypothesis tests about the world.

Real-world applications include forecasting sales from marketing spend, estimating housing prices from square footage and neighbourhood features, pricing insurance policies from customer attributes, modelling how drug dosage affects blood glucose, projecting energy demand from temperature, and, more generally, any situation where a continuous outcome has to be explained by a mix of numerical and categorical drivers.

2. Business Problem: Predicting Medical Insurance Costs

Imagine you work as a data scientist for a mid-sized health insurance company. The underwriting team currently sets premiums using a rule book that was designed a decade ago: they look up the applicant's age band, apply a smoker surcharge, add a regional multiplier, and quote a number. It works, but it is coarse. Two 40-year-old non-smokers with very different body-mass-index profiles pay the same premium, even though one is dramatically more likely to file a large claim. Some customers are being overcharged and defect to competitors; others are being undercharged and quietly bleed the company money.

The business objective is to estimate a customer's expected annual medical charges before a policy is issued, using the small handful of attributes captured during application. A more accurate estimate lets the pricing team set fair, competitive premiums, tightens the loss ratio, and gives the marketing team a defensible story about why prices differ. Success will be measured with three metrics: root mean squared error on a held-out set, the R-squared of the model, and — just as important — the interpretability of the coefficients, because the regulator has to be able to read the model and understand it.

OLS is the right first tool for this problem. The target — annual charges — is a continuous number. The predictors are a modest mix of numerical and categorical variables. The stakeholders need a model they can explain, not a black box. And the training set is small enough (a few thousand rows) that a linear model will generalise more reliably than a deep neural network.

3. Understanding the Dataset

We will use the well-known Medical Cost Personal Dataset, often called the Insurance dataset. It contains 1,338 rows and seven columns. Each row is one policyholder. Six columns are predictors; the seventh, charges, is the target we want to model.

• age — the primary policyholder's age in years. Numeric. Older customers typically incur higher medical costs.
• sex — biological sex, either male or female. Categorical. Some medical conditions have sex-linked risk profiles.
• bmi — body mass index, weight in kilograms divided by height in metres squared. Numeric. A well-documented driver of chronic-disease risk.
• children — number of dependents covered by the policy. Integer count. Affects the size of the covered pool.
• smoker — yes or no. Categorical. The single largest driver of medical charges in this dataset.
• region — one of northeast, northwest, southeast, southwest. Categorical. Absorbs cost-of-care differences across geographies.
• charges — annual medical expenses billed to the insurer, in US dollars. Continuous. This is our target variable.

A quick preview of the first three rows would look roughly like: (19, female, 27.9, 0, yes, southwest, 16884.92), (18, male, 33.77, 1, no, southeast, 1725.55), (28, male, 33.0, 3, no, southeast, 4449.46). Even from three rows the smoker column looks suspicious — the smoking teenager is paying ten times what her non-smoking peers pay. That intuition will show up loudly in the model.

4. Exploratory Data Analysis

Before we fit anything, we look. EDA is not an optional appetiser — it is where you learn the personality of your data, catch broken columns, and form the hypotheses that the model will later confirm or refute.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Download the dataset from a public HTTP URL, cache it locally, then load it.
import os
CSV_SOURCE = "https://raw.githubusercontent.com/stedy/Machine-Learning-with-R-datasets/master/insurance.csv"
CSV_URL = "insurance.csv"
if not os.path.exists(CSV_URL):
    try:
        # In the browser (Pyodide), stdlib urllib can't do https — use pyfetch.
        from pyodide.http import pyfetch
        resp = await pyfetch(CSV_SOURCE)
        with open(CSV_URL, "wb") as f:
            f.write(await resp.bytes())
    except ImportError:
        import urllib.request
        urllib.request.urlretrieve(CSV_SOURCE, CSV_URL)
df = pd.read_csv(CSV_URL)

# 1. Structural checks
print(df.shape)              # (1338, 7)
print(df.dtypes)             # dtypes per column
print(df.isnull().sum())     # missing values per column
print(df.duplicated().sum()) # exact duplicate rows

# 2. Summary statistics
print(df.describe(include="all").T)

# 3. Target distribution
plt.figure(figsize=(8, 4))
sns.histplot(df["charges"], bins=40, kde=True)
plt.title("Distribution of Annual Medical Charges")
plt.xlabel("Charges (USD)"); plt.ylabel("Count")
plt.show()

# 4. Categorical breakdown of the target
plt.figure(figsize=(8, 4))
sns.boxplot(x="smoker", y="charges", data=df)
plt.title("Charges by Smoking Status")
plt.show()

# 5. Numeric relationship with the target
sns.scatterplot(x="bmi", y="charges", hue="smoker", data=df)
plt.title("Charges vs BMI, coloured by smoking status")
plt.show()

# 6. Correlation among numeric columns
num_cols = ["age", "bmi", "children", "charges"]
plt.figure(figsize=(6, 5))
sns.heatmap(df[num_cols].corr(), annot=True, cmap="coolwarm", fmt=".2f")
plt.title("Correlation Heatmap")
plt.show()

Reading these plots in order tells a coherent story. The histogram of charges is heavily right-skewed — most policyholders cost the insurer a few thousand dollars a year, but a long tail runs above forty thousand. That skew matters: OLS does not require the target itself to be normal, but extreme tails can distort residuals and inflate error metrics. The box plot of charges by smoker is the single most informative chart in the entire notebook. The median non-smoker sits near eight thousand dollars; the median smoker sits above thirty thousand. That gap of roughly twenty-three thousand dollars is the smoker surcharge in raw form.

The scatter of charges against BMI, coloured by smoking status, reveals something even more interesting: the two groups behave differently. For non-smokers, BMI barely nudges charges — the cloud is nearly flat. For smokers, charges climb sharply with BMI, with a visible kink around a BMI of thirty (the clinical obesity threshold). This is called an interaction, and it is a hint that a plain additive linear model may be leaving accuracy on the table. The correlation heatmap on the remaining numerics is unsurprising: age correlates moderately with charges (around 0.30), BMI weakly (around 0.20), and children hardly at all.

Business insight from EDA: the pricing team's intuition that smoking dominates is correct, but there is a second-order effect worth surfacing — obese smokers are a distinctly high-risk cohort and deserve their own actuarial treatment. Common mistake to avoid: do not delete the high-charge outliers. In an insurance context, those tails are not noise — they are the exact events the premium is meant to cover.

5. Data Preprocessing

OLS is comparatively forgiving about scale — a linear model with an intercept absorbs constant shifts and unit changes without any loss of accuracy. What it will not tolerate is text data. Every categorical column must be turned into numbers, and every missing value must be either imputed or removed. In our dataset there are no missing values and no duplicates, so the preprocessing focuses on encoding.

from sklearn.model_selection import train_test_split

# One-hot encode the categorical columns.
# drop_first=True avoids the dummy-variable trap: for k categories we
# create k-1 dummies so the design matrix stays full-rank.
df_enc = pd.get_dummies(
    df,
    columns=["sex", "smoker", "region"],
    drop_first=True,
    dtype=int,
)

# Split predictors (X) and target (y).
X = df_enc.drop(columns=["charges"])
y = df_enc["charges"]

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.20, random_state=42
)
print(X_train.shape, X_test.shape)

A few decisions deserve comment. We used one-hot encoding rather than label encoding because the categorical variables have no natural order — mapping northeast to 0, northwest to 1 and so on would inject a fake ranking into the model. We passed drop_first=True to sidestep perfect multicollinearity between the dummy columns and the intercept. We did not scale the numeric columns; OLS coefficients simply pick up the appropriate units. And we set a fixed random_state so the train/test split is reproducible for every teammate who re-runs the notebook.

Common mistakes here include leaking test information into training (for example, computing summary statistics on the full dataset before the split), forgetting to drop one dummy per categorical variable, and dropping a categorical column altogether because it looked messy. Each of these silently damages either the estimated coefficients or the honest evaluation of the model.

6. Mathematical Foundation of OLS

Linear regression assumes that the target y can be written as a weighted sum of the predictors, plus an intercept, plus a random error term. For a single observation with p predictors:

Here beta zero is the intercept — the value of y when every predictor is zero — and each beta j is the slope for predictor j: the amount y changes when x j increases by one unit, holding every other predictor fixed. The epsilon term captures everything the linear model cannot: measurement noise, omitted variables, genuinely random variation.

Stack all n observations together and the equation becomes matrix arithmetic. Let X be an n-by-(p+1) matrix whose first column is a column of ones (for the intercept) and whose remaining columns hold the predictor values. Let y be a column vector of length n, and beta a column vector of length p+1. Then the model is simply:

For any candidate beta, the model's prediction on the training data is X beta and the residual — the error the model makes on row i — is r_i = y_i minus the i-th prediction. OLS chooses the beta that minimises the sum of squared residuals, which we call the residual sum of squares:

Why squared error and not absolute error? Three reasons. First, squaring is differentiable everywhere, which makes the optimisation clean. Second, squaring penalises big errors much more than small ones, which matches the intuition that being off by ten thousand dollars is more than twice as bad as being off by five thousand. Third — and this is deep — if the errors are normally distributed, the least-squares estimator coincides with the maximum-likelihood estimator, so we are secretly doing statistics of the most principled kind.

To find the minimum, take the gradient of RSS with respect to beta, set it to zero, and solve. The algebra takes about four lines and produces one of the most elegant results in applied mathematics — the Normal Equation:

Geometrically, X beta lives in the column space of X — the set of all vectors we can build by combining the predictor columns. The best beta is the one that makes X beta the orthogonal projection of y onto that column space. Statistically, under the Gauss–Markov assumptions we discuss next, this estimator is BLUE — the Best Linear Unbiased Estimator — meaning it has the smallest variance among all linear unbiased estimators of beta. That is a remarkable optimality guarantee for such a simple formula.

7. The Five OLS Assumptions

OLS is optimal when five assumptions hold. Real data almost never satisfies all five perfectly, so the practical skill is knowing which violations you can live with, which distort your estimates, and which invalidate your inference.

Linearity

The relationship between predictors and target must be linear in the coefficients. Check with a residual-versus-fitted plot: if you see a systematic curve, the linearity assumption is broken. Remedies include adding polynomial terms, log-transforming skewed predictors, or switching to a genuinely non-linear model such as a gradient-boosted tree.

Independence of Errors

Each residual must be independent of every other residual. In cross-sectional data (like our insurance dataset) this usually holds. In time-series data it almost never does — today's error tends to look like yesterday's error. The Durbin–Watson statistic near two signals independence; values much below two indicate positive autocorrelation. Remedy: add lag terms or move to a time-aware model such as ARIMA.

Homoscedasticity

The variance of the residuals must be constant across the range of fitted values. On a residual-versus-fitted plot this looks like a horizontal band; a fan shape signals heteroscedasticity. Formal tests include Breusch–Pagan and White. When violated, coefficient estimates remain unbiased but standard errors are wrong, so p-values and confidence intervals become unreliable. Remedies: log-transform the target, use weighted least squares, or use heteroscedasticity-consistent (HC3) standard errors.

Normality of Residuals

The residuals — not the predictors, not the target — should be approximately normally distributed. This assumption powers the t-tests and F-tests on the coefficients. Check with a Q–Q plot or the Shapiro–Wilk test. With large samples (say n above a few hundred) the central limit theorem is forgiving and mild deviations from normality rarely change conclusions.

No Multicollinearity

Predictors should not be highly linearly dependent on each other. If two columns carry essentially the same information, the matrix X-transpose-X becomes near-singular, the inverse blows up, and the coefficient estimates become wildly unstable — flipping sign when you add or remove a single row. The standard diagnostic is the Variance Inflation Factor. A VIF above five is a warning; above ten is a serious problem. Remedies include dropping one of the offending predictors, combining them into a single index, or switching to Ridge regression, which is essentially OLS with a built-in regulariser for exactly this problem.

import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.diagnostic import het_breuschpagan
from scipy import stats

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

# Multicollinearity
vif = pd.DataFrame({
    "feature": X_train_c.columns,
    "VIF": [variance_inflation_factor(X_train_c.values, i)
            for i in range(X_train_c.shape[1])],
})
print(vif)

# Homoscedasticity (Breusch-Pagan)
bp = het_breuschpagan(model.resid, model.model.exog)
print("Breusch-Pagan p-value:", bp[1])

# Normality of residuals (Shapiro-Wilk)
print("Shapiro-Wilk p-value:", stats.shapiro(model.resid).pvalue)

8. Building the OLS Model

The scikit-learn LinearRegression class is fine for prediction, but statsmodels is the right choice when you also want inference — standard errors, p-values, confidence intervals. The API is deliberately close to R's lm() function.

import statsmodels.api as sm

# statsmodels does not add an intercept automatically.
X_train_const = sm.add_constant(X_train)
X_test_const  = sm.add_constant(X_test)

model = sm.OLS(y_train, X_train_const).fit()
print(model.summary())

The summary printed by statsmodels is dense but every number earns its place. Read it in three passes: overall fit, individual coefficients, then diagnostics.

Overall fit. R-squared reports the fraction of variance in charges the model explains — typically around 0.75 on this dataset, meaning three-quarters of the variation is captured. Adjusted R-squared corrects R-squared for the number of predictors, penalising you for adding variables that do not carry their weight; if it barely differs from R-squared, your predictors are pulling their weight. The F-statistic and its p-value ask a single yes/no question: is at least one predictor useful? A p-value below 0.05 says yes.

Individual coefficients. For each row, the coefficient is the estimated slope, the standard error is our uncertainty about that slope, the t-statistic is coefficient divided by standard error, and the p-value is the probability of seeing a t-statistic that extreme if the true coefficient were zero. The 95 percent confidence interval is the range of values consistent with the data. In practical terms, a smoker coefficient around 23,800 with a p-value under 0.001 and a confidence interval that comfortably excludes zero says: after controlling for every other variable, being a smoker adds roughly twenty-four thousand dollars to expected annual charges, and we are very confident that this is not luck.

Model-selection diagnostics. AIC (Akaike Information Criterion) and BIC (Bayesian Information Criterion) are useful for comparing models: lower is better. Both trade goodness of fit against complexity, with BIC penalising complexity more aggressively. The residual standard error tells you the typical size of a prediction error in the units of the target — around six thousand dollars for this model. Degrees of freedom are just n minus (p+1) — the amount of information left after fitting the coefficients.

9. Model Diagnostics

The regression summary tells you whether the arithmetic worked. The diagnostic plots tell you whether the model is reasonable. Skipping this step is the single most common mistake in applied regression.

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm

fitted   = model.fittedvalues
resid    = model.resid
influence = model.get_influence()
std_resid = influence.resid_studentized_internal
leverage  = influence.hat_matrix_diag
cooks     = influence.cooks_distance[0]

fig, axes = plt.subplots(2, 2, figsize=(12, 9))

# Residuals vs Fitted — checks linearity & homoscedasticity
axes[0, 0].scatter(fitted, resid, alpha=0.4)
axes[0, 0].axhline(0, color="red", lw=1)
axes[0, 0].set_title("Residuals vs Fitted")

# Q-Q plot — checks normality
sm.qqplot(resid, line="45", fit=True, ax=axes[0, 1])
axes[0, 1].set_title("Normal Q-Q")

# Scale-Location — checks homoscedasticity
axes[1, 0].scatter(fitted, np.sqrt(np.abs(std_resid)), alpha=0.4)
axes[1, 0].set_title("Scale-Location")

# Cook's Distance — flags influential points
axes[1, 1].stem(np.arange(len(cooks)), cooks, markerfmt=",")
axes[1, 1].set_title("Cook's Distance")

plt.tight_layout(); plt.show()

Reading the plots: a healthy residuals-versus-fitted plot looks like a shapeless cloud around the zero line. On this insurance dataset you will actually see two distinct bands of points — that is the smoker/non-smoker split showing up in the residuals, and it is the strongest hint that adding a smoker-by-BMI interaction term will materially improve the model. A well-behaved Q-Q plot has the points hugging the 45-degree line; deviations at the tails indicate that extreme charges are heavier-tailed than a normal distribution predicts. The scale-location plot should be a horizontal band; a rising trend signals heteroscedasticity. Cook's distance flags individual rows whose removal would meaningfully shift the coefficients — any point above 4/n deserves a look, but do not delete outliers reflexively.

10. Model Evaluation

Evaluation happens on the held-out test set, not the training set. The five standard regression metrics answer complementary questions.

MAE is the average absolute error, in the units of the target, and is robust to outliers. MSE squares the errors, so it punishes big misses harder — the natural loss for OLS itself. RMSE is the square root of MSE, back in target units, and is the metric non-technical stakeholders usually understand fastest. R-squared measures how much of the variance the model explains; a value of one is perfect, zero means the model does no better than predicting the mean. Adjusted R-squared subtracts a penalty for the number of predictors, discouraging you from cramming useless features into the model.

from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
import numpy as np, matplotlib.pyplot as plt

y_pred = model.predict(X_test_const)

mae  = mean_absolute_error(y_test, y_pred)
mse  = mean_squared_error(y_test, y_pred)
rmse = np.sqrt(mse)
r2   = r2_score(y_test, y_pred)
n, p = X_test.shape
adj_r2 = 1 - (1 - r2) * (n - 1) / (n - p - 1)

print(f"MAE : {mae:,.2f}")
print(f"RMSE: {rmse:,.2f}")
print(f"R^2 : {r2:.4f}   Adj R^2: {adj_r2:.4f}")

# Actual vs predicted
plt.figure(figsize=(7, 6))
plt.scatter(y_test, y_pred, alpha=0.5)
lims = [y_test.min(), y_test.max()]
plt.plot(lims, lims, "r--")
plt.xlabel("Actual charges"); plt.ylabel("Predicted charges")
plt.title("Actual vs Predicted")
plt.show()

On this dataset you should see RMSE around six thousand dollars and R-squared near 0.78 on the test set. The actual-versus-predicted scatter should hug the diagonal for typical customers, with two visible groups peeling upward — smokers and, especially, obese smokers whose costs the plain additive model systematically under-predicts. That systematic error is exactly the story the diagnostics warned us about.

11. Business Interpretation

This is where the analyst earns her salary. A statistical result is only valuable when it changes a decision. Reading the coefficients back in plain business language on the insurance model:

• Age contributes roughly $250 in additional annual charges per year of age, holding everything else constant. A 60-year-old customer is therefore expected to cost about $10,000 more per year than a 20-year-old with the same profile.
• BMI adds around $340 per unit. Moving from a healthy BMI of 22 to an obese BMI of 32 adds roughly $3,400 to expected annual charges.
• Being a smoker adds approximately $23,800 in expected annual charges — by an enormous margin the single most important variable in the model, and completely consistent with the actuarial literature.
• Each additional dependent adds a small positive amount (around $400) but the coefficient is statistically weak — this variable barely earns its place.
• Region coefficients are small (a few hundred dollars) but statistically significant for the southeast, capturing regional cost-of-care differences.
• Sex is not statistically significant once age, BMI, smoking status and region are controlled for. That is a genuinely useful finding: it tells the underwriting team they can drop sex as a pricing variable without harming accuracy.

Recommended action for the pricing team: keep the smoker surcharge as the dominant lever, add an obese-smoker cross term for actuarial fairness, treat age as a smooth continuous factor rather than crude age bands, and drop sex from the rating factors. The regulator will appreciate a model whose coefficients each tell a defensible story.

12. Complete End-to-End Project

The following script stitches every step together into a single reproducible pipeline you can run start-to-finish on the insurance dataset.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
import joblib
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

# 1. Download the dataset from a public HTTP URL, cache it locally, then load it.
import os
CSV_SOURCE = "https://raw.githubusercontent.com/stedy/Machine-Learning-with-R-datasets/master/insurance.csv"
CSV_URL = "insurance.csv"
if not os.path.exists(CSV_URL):
    try:
        # In the browser (Pyodide), stdlib urllib can't do https — use pyfetch.
        from pyodide.http import pyfetch
        resp = await pyfetch(CSV_SOURCE)
        with open(CSV_URL, "wb") as f:
            f.write(await resp.bytes())
    except ImportError:
        import urllib.request
        urllib.request.urlretrieve(CSV_SOURCE, CSV_URL)
df = pd.read_csv(CSV_URL)

# 2. Quick EDA
print(df.shape, df.isnull().sum().sum(), "duplicates:", df.duplicated().sum())

# 3. Feature engineering: interaction term for obese smokers
df["bmi30"] = (df["bmi"] >= 30).astype(int)
df_enc = pd.get_dummies(
    df, columns=["sex", "smoker", "region"], drop_first=True, dtype=int
)
df_enc["smoker_bmi30"] = df_enc["smoker_yes"] * df_enc["bmi30"]

# 4. Split
X = df_enc.drop(columns=["charges"])
y = df_enc["charges"]
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.20, random_state=42
)

# 5. Fit OLS with statsmodels
X_train_c = sm.add_constant(X_train)
X_test_c  = sm.add_constant(X_test)
model = sm.OLS(y_train, X_train_c).fit()
print(model.summary())

# 6. Evaluate
pred = model.predict(X_test_c)
print(f"MAE : {mean_absolute_error(y_test, pred):,.2f}")
print(f"RMSE: {np.sqrt(mean_squared_error(y_test, pred)):,.2f}")
print(f"R^2 : {r2_score(y_test, pred):.4f}")

# 7. Predict on a new customer
new_customer = pd.DataFrame([{
    "age": 35, "bmi": 31.2, "children": 2, "bmi30": 1,
    "sex_male": 1, "smoker_yes": 1,
    "region_northwest": 0, "region_southeast": 1, "region_southwest": 0,
    "smoker_bmi30": 1,
}])
new_customer_c = sm.add_constant(new_customer, has_constant="add")
print("Predicted annual charge:",
      f"$ {model.predict(new_customer_c).iloc[0]:,.2f}")

# 8. Persist and reload
joblib.dump(model, "ols_insurance.joblib")
reloaded = joblib.load("ols_insurance.joblib")

Adding the obese-smoker interaction typically lifts R-squared from around 0.75 to above 0.85, cuts RMSE by roughly a third, and makes the actual-versus-predicted plot noticeably tighter. That is what a well-diagnosed model looks like: the residual plot did not just complain — it pointed straight at the fix.

13. Common Mistakes

• Reading correlation as causation. A large coefficient means the two variables move together after controlling for the others in the model — nothing about cause and effect. That takes an experiment or a well-designed quasi-experimental study.
• Blindly trusting R-squared. A model can hit R-squared of 0.99 by memorising noise on a small dataset. Always evaluate on held-out data and inspect residuals.
• Ignoring multicollinearity. Two near-identical predictors will produce two enormous, opposite coefficients that cancel out. The predictions still work; the interpretation is nonsense.
• Interpreting a p-value as the probability the null hypothesis is true. It is not. It is the probability of the observed data (or something more extreme) assuming the null is true.
• Data leakage: fitting scalers, encoders or imputers on the full dataset before splitting. This subtly inflates every metric you compute.
• Deleting outliers because they are inconvenient. In insurance, healthcare and finance, the outliers are frequently the point of the exercise.
• Failing to check assumptions and then quoting coefficient p-values as if they were gospel. Standard errors are only trustworthy when the assumptions hold.
• Overfitting by throwing every column and every interaction into the model. Use adjusted R-squared, AIC/BIC and cross-validation to control complexity.

14. Interview Questions and Answers

Foundations

• Q1. What does OLS actually optimise? The sum of squared residuals between observed and predicted values.
• Q2. Why squared error rather than absolute error? Differentiability, harsher penalty on large errors, and equivalence with maximum likelihood under Gaussian noise.
• Q3. State the Normal Equation. beta-hat equals (X-transpose X) inverse times X-transpose y.
• Q4. What is the Gauss–Markov theorem? Under linearity, independence, homoscedasticity and zero-mean errors, OLS is the Best Linear Unbiased Estimator of beta — smallest variance among all linear unbiased estimators.
• Q5. Difference between OLS and Gradient Descent? OLS is a closed-form solution; gradient descent is an iterative optimiser that reaches the same solution and is preferred when X-transpose X is huge or singular.
• Q6. Why must we add a column of ones to X? To let the model fit a non-zero intercept.
• Q7. When would the Normal Equation fail? When X-transpose X is singular — caused by perfect multicollinearity or by having more predictors than observations.

Assumptions and Diagnostics

• Q8. Name the five OLS assumptions. Linearity, independence of errors, homoscedasticity, normality of residuals, and no multicollinearity.
• Q9. What is heteroscedasticity and how do you detect it? Non-constant residual variance; detect visually with residuals-vs-fitted or scale-location plots, formally with Breusch–Pagan or White tests.
• Q10. Fix for heteroscedasticity? Log-transform the target, use weighted least squares, or report HC3 robust standard errors.
• Q11. What is multicollinearity? Strong linear dependence between predictors. Detect with the Variance Inflation Factor; VIF above 10 is a red flag.
• Q12. What is VIF? The factor by which the variance of a coefficient is inflated because that predictor is correlated with the others. VIF equals 1 divided by (1 minus R-squared of regressing that predictor on the rest).
• Q13. What does the Durbin–Watson statistic measure? Autocorrelation in residuals; values near 2 indicate no autocorrelation.
• Q14. How do you check normality of residuals? Q-Q plot, Shapiro–Wilk test, Anderson–Darling test.
• Q15. What is Cook's distance? A measure of how much the fitted values would change if a single observation were removed — flags influential points.

Interpretation and Evaluation

• Q16. Interpret a coefficient of 250 on age. Holding all other predictors constant, a one-year increase in age is associated with a $250 increase in expected charges.
• Q17. What is the difference between R-squared and adjusted R-squared? Adjusted R-squared penalises R-squared for the number of predictors, so it does not automatically rise when useless features are added.
• Q18. Why can a high R-squared still indicate a poor model? Because R-squared measures variance explained on the training data, not predictive accuracy on new data, and it does not detect assumption violations.
• Q19. What is the F-statistic in the regression summary? A joint test of whether at least one predictor coefficient is non-zero.
• Q20. Difference between AIC and BIC? Both trade fit against complexity; BIC penalises complexity more heavily and is preferred when the true model is believed to be sparse.
• Q21. What does a p-value of 0.03 on a coefficient mean? If the true coefficient were zero, we would observe an estimate this extreme only 3% of the time.
• Q22. When is MAE preferable to RMSE? When large errors should not be over-penalised, or when the target distribution has heavy outliers.

Practical and Scenario Based

• Q23. When would you use OLS instead of Ridge? When you have low multicollinearity, plenty of data relative to predictors, and inference matters more than regularised prediction.
• Q24. When would Ridge or Lasso beat OLS? With many correlated predictors, or when you need automatic feature selection (Lasso).
• Q25. Your R-squared is 0.99 on train but 0.30 on test. What is happening? Textbook overfitting — reduce complexity, add regularisation, cross-validate.
• Q26. Your coefficient on a feature flipped sign when you added another feature. Why? Almost certainly multicollinearity between the two features.
• Q27. Residuals-vs-fitted plot shows a U-shape. What now? Non-linearity — add polynomial terms, log-transform, or use a non-linear model.
• Q28. How do you handle categorical features with 100 levels? Target encoding, frequency encoding or embedding rather than one-hot, to keep the design matrix tractable.
• Q29. Can OLS be used for classification? Directly, no — probabilities are not bounded to [0,1]. Use logistic regression, which is OLS's classification cousin.
• Q30. Explain the bias-variance trade-off in the OLS context. Adding features reduces bias but increases variance; regularisation trades a little bias for a lot of variance reduction.

15. Practice Exercises

• Refit the model after removing observations whose charges exceed $50,000. How much do the coefficients move? Does R-squared genuinely improve or just look better?
• Add polynomial terms in BMI (BMI squared, BMI cubed). Does adjusted R-squared improve, or are the extra terms just noise?
• Fit Ridge and Lasso on the same design matrix using sklearn. Compare coefficients, RMSE and which features Lasso zeros out.
• Perform backward elimination: drop the least significant predictor one at a time until every remaining p-value is below 0.05. Compare the reduced model to the full one on the test set.
• Try a 60/40, 70/30 and 80/20 train/test split with the same random seed. How much does test RMSE bounce around? What does that tell you about the reliability of a single split?
• Predict expected annual charges for three fictional customers of your own design and defend the numbers to a hypothetical underwriting manager.

16. Conclusion

Ordinary Least Squares is old, unglamorous, and often the best tool for the job. In one fitted object it hands you a prediction machine, a set of hypothesis tests about the world, and a diagnostic apparatus rich enough to tell you exactly how the model is failing. On the insurance dataset it identified smoking as the dominant cost driver, quantified the age and BMI effects in defensible dollars, quietly retired sex as a rating factor, and — through the diagnostic plots — pointed us toward the interaction term that lifted R-squared into the mid-eighties.

The workflow generalises. Load and inspect the data. Encode and split. Fit an OLS model in statsmodels. Read the summary from top to bottom. Check the five assumptions with residual plots, Q-Q plots, VIF and Breusch–Pagan. Evaluate on a held-out set with MAE, RMSE and R-squared. Translate the coefficients into business language, and be honest about what the model cannot see. When those checks fail — when residuals fan out, when a coefficient flips sign, when a diagnostic plot curls in a suspicious way — the failure is telling you something specific about the world your data came from, and the fix is almost always to enrich the model rather than abandon it.

Once OLS feels natural, the doors it opens are enormous. Multiple linear regression is already what you have been doing. Polynomial regression is a two-line extension. Ridge, Lasso and Elastic Net are OLS with a regularisation term bolted on. Logistic regression swaps the identity link for the logit and inherits almost the same summary table. Generalised linear models extend the idea to Poisson counts, gamma-distributed costs and beyond. Learn OLS well and the rest of applied statistics stops feeling like a zoo of unrelated methods and starts looking like variations on a single, beautiful theme.