6.1 OLS, Coefficients, and the Design Matrix: Exactly What Is a Linear Model Estimating?

Sample differences, interval estimates, and test statistics provide different kinds of evidence; one p-value cannot stamp a conclusion “reliable.” To describe how several continuous factors vary with an outcome, place them in one regression design matrix.

“If resources increase by 10%, how much will completion improve?” A regression can describe conditional associations, but the effect of an intervention requires a suitable design or identifying assumptions. Decide whether the task is prediction, association or causal estimation before choosing the columns.

This lesson's objectives

• Write a multiple linear regression model using a design matrix;
• Explain the OLS objective and the conditional interpretation of the coefficients;
• Properly handle categories, interactions, and nonlinear terms;
• Distinguish predictive models, statistical inference, and causal questions.

1. Linearity means linearity with respect to parameters

Multiple linear regression:

$$Y_i={\beta }_0+{\beta }_1{X}_{i1}+\cdots +{\beta }_p{X}_{ip}+{\epsilon }_i.$$

Matrix form:

$$\mathbf{y}=X\mathit{\beta }+\mathit{\epsilon }.$$

"Linear" means linearity in the coefficients $\beta$. The model can include $x^2$, $\log x$, splines, and interactions like $xz$, as long as the coefficients multiply fixed transformed columns and enter linearly.

2. OLS minimizes the sum of squared residuals

$$\hat{\mathit{\beta }}=\arg \underset{\mathit{\beta }}{min}\Vert \mathbf{y}-X\mathit{\beta }{\Vert }_2^2.$$

Under full rank, the normal equations can be written, but numerical implementations typically use QR or SVD, avoiding the explicit computation of $(X^{T}X{)}^{-1}$. Volume 10 already explained how forming $X^{T}X$ worsens the condition number.

The columns must have full rank for a unique coefficient vector. A rank-deficient design can still have unique fitted values on the training rows, while individual coefficients remain unidentified; a solver returning numbers does not settle their interpretation.

Squared loss gives higher weight to large residuals, making OLS sensitive to influential points. It does not assume that $Y$ itself must be normal; assumptions about the error distribution primarily affect finite-sample inference.

3. Coefficients are the associations "when other columns are fixed"

Model:

$$duration={\beta }_0+{\beta }_{1}team\mathrm{\_}size+{\beta }_{2}difficulty+\epsilon .$$

${\beta }_1$ represents the conditional difference in average duration when team size increases by 1, holding all other variables like difficulty constant.

It doesn't automatically imply that adding a person to the same task causes this change. Insufficient measurement of task difficulty, resource allocation mechanisms, and reverse causality could all introduce confounding factors. Causal interpretation requires either an appropriate study design or additional identifying assumptions.

4. Interpret the Intercept Within Its Meaningful Range

${\beta }_0$ is the prediction when all numeric features are 0 and the category is in the reference group. If team size = 0 is not in the data range, the intercept is merely a parameterization necessity and lacks real-world meaning.

Centering continuous variables puts the intercept in line with the mean:

The synthetic fixture below runs with pandas, NumPy and statsmodels 0.15. Run the Python blocks in order. Both numeric predictors are centered using training means, so the reference-category intercept is evaluated at those numeric means.

import numpy as np
import pandas as pd

rng = np.random.default_rng(42)
training = pd.DataFrame({
    "team_size": rng.integers(2, 9, 48),
    "resources_used": rng.uniform(10, 50, 48),
    "mission_type": np.tile(["survey", "rescue", "remote"], 16),
})
training["duration_minutes"] = (
    30 - 2 * training["team_size"] + 0.4 * training["resources_used"]
    + training["mission_type"].map({"survey": 0, "rescue": 2, "remote": 6})
    - 0.5 * training["team_size"] * training["mission_type"].eq("rescue")
    + rng.normal(0, 2, len(training))
)
centers = training[["team_size", "resources_used"]].mean()
for name in centers.index:
    training[name + "_centered"] = training[name] - centers[name]

For evaluation or online prediction, the mean can only be estimated during training and then stored.

5. Categorical variables require a reference group or constraint

Nominal task types should not be coded as 1,2,3 and treated as equally spaced numeric values. For an unregularized OLS model with an intercept, use one omitted reference category or another explicit constraint to remove the full dummy set’s linear dependence.

Intercept: reference-type prediction when all other numeric columns are zero
type_rescue: differences in conditions between rescue and reference
type_remote: the conditional difference between remote and reference

Changing the reference group doesn't affect the fitted values, but does change the interpretation of individual coefficients. When reporting, list the reference group and not just show the column names.

6. Interactions Let Slopes Vary with Conditions

duration ~ resources + difficulty + resources:difficulty

Model:

$$Y={\beta }_0+{\beta }_{1}R+{\beta }_{2}D+{\beta }_{3}RD+\epsilon .$$

The slope with respect to resources is ${\beta }_1+{\beta }_{3}D$, no longer fixed at ${\beta }_1$. Interpreting main effects alone when an interaction is included would be incorrect.

The hierarchy principle typically recommends retaining main effects when interactions are present, even if individual p-values are not significant; model structure should arise from the research question, not from stepwise removal until all terms are significant.

7. Use the formula interface to clearly define the model

import statsmodels.formula.api as smf

formula = (
    'duration_minutes ~ resources_used_centered + team_size_centered'
    ' * C(mission_type, Treatment(reference="survey"))'
)
model = smf.ols(formula, data=training, missing="raise").fit()
if np.linalg.matrix_rank(model.model.exog) != model.model.exog.shape[1]:
    raise ValueError("The design does not identify every coefficient")
print(model.params)
print(model.nobs, model.model.exog_names)

The formula interface makes it easy to handle factors and interactions, but specific coding, reference levels, missing value removal, and robust covariance options need to be checked. summary() is not an automatic correctness guarantee.

Here * expands to both main effects and their interaction. The coefficient for a task category is its difference from survey at the centered team-size value zero, not at every team size. Inspect the design-column names and row count; unseen categories at prediction time also need an explicit policy.

8. $R^2$ Measures the relative size of sample squared error compared to a baseline

$$R^2=1-\frac{\sum _i(y_i-{\hat{y}}_i{)}^2}{\sum _i(y_i-\overline{y}{)}^2}.$$

For nested, unregularized OLS models fitted to the same rows and target with an intercept, adding columns cannot increase the minimized training SSE, so ordinary training $R^2$ cannot decrease. Changing the rows or adding regularization removes that guarantee. In the usual test-set formula, $\overline{y}$ is the evaluation set’s mean. A negative test $R^2$ means worse squared error than that mean, not necessarily worse than a deployable training-mean baseline.

from sklearn.metrics import r2_score

y_test = np.array([10., 12.])
predicted = np.array([10., 10.])
training_mean = 0.0
print(r2_score(y_test, predicted))  # -1.0
print(np.sum((y_test - predicted) ** 2))  # 4.0
print(np.sum((y_test - training_mean) ** 2))  # 244.0

This deliberately shifted example has negative $R^2$ while outperforming the training-mean baseline. The evaluation mean is a scoring reference, not a prediction available in advance. Constant evaluation targets make the raw ratio undefined; scikit-learn R² documents its finite-value replacement policy.

$R^2$ does not imply unbiased coefficients, valid residuals, causal interpretation, or sufficiently small prediction errors. MAE/RMSE and baseline should be reported on the target scale.

9. Result types might not be suitable for OLS

If each row's outcome is a binary success/failure, an OLS linear probability model may predict values less than 0 or greater than 1 and exhibit structural heteroskedasticity; logistic regression typically fits better within the probability bounds.

If the goal is to aggregate success_rate, different denominators lead to different variances, requiring weights or a binomial model. Count data, positive long tails, and censored time points each have more appropriate model families.

OLS can still serve as a simple baseline or be used for proportional outcomes under clearly defined conditions, but diagnostic evaluation and acknowledgment of limitations are required, defaulting to OLS merely because it's familiar with the API is not acceptable.

Common Misconceptions

• A coefficient represents a variable's causal impact on the outcome: Observed regression defaults to conditional association only.
• Linear models cannot capture curved relationships: basis functions can be added, keeping the parameters linear.
• A high $R^2$ indicates a correct model: leakage and misspecification can still produce high values.
• Outcomes between 0 and 1 can be regressed without further checks: the denominator and distribution model need to be handled.

Practice

1. Write an "all other columns fixed" interpretation for each coefficient in the task duration model.
2. Change the reference group in the category, verify that the fitted values remain unchanged but the coefficients do change.
3. Introduce a resource-difficulty interaction and derive the resource slope for different difficulty levels.
4. Construct a model with high $R^2$ caused by time leakage.

Summary

OLS estimates the conditional mean structure by minimizing squared residuals. The design matrix determines the interpretation of coefficients, while categories, interactions, and transformations dictate what the model can express; prediction, association, and causation are three distinct objectives.

Next lesson: diagnosing where the model goes wrong, nonlinearity, heteroskedasticity, correlated errors, high-leverage points, and multicollinearity each affect predictions and inference in distinct ways.