Multiple Regression: An Introduction

Simple regression fits a line, $\hat{Y} = b_0 + b_1 X$, describing how $Y$ changes with a single predictor $X$. Multiple regression generalizes this to a plane (with two predictors) or a hyperplane (with many):

Here $Y$ is the response variable, and $X_1, \ldots, X_k$ are explanatory variables that may include continuous measurements, transformed variables, or indicator (dummy) variables for categorical groups. Each partial regression coefficient $b_j$ estimates the expected change in $Y$ for a one-unit increase in $X_j$, holding all other $X$ variables in the model fixed. This qualifier — "holding constant" — is the entire point of the method. It is what separates a partial coefficient from the crude, marginal association you'd get by regressing $Y$ on $X_j$ alone. The coefficients are still chosen by ordinary least squares (OLS): the values of $b_0, b_1, \ldots, b_k$ that minimize the sum of squared residuals,

across all $n$ observations simultaneously. Geometrically, with two predictors this fits the best-fitting plane through a cloud of points in three-dimensional space (one axis for $Y$, one for $X_1$, one for $X_2$); with $k$ predictors it fits a hyperplane in $(k+1)$-dimensional space. The residual for each point is now the vertical distance to this hyperplane, not to a line. A crucial conceptual bridge: a partial coefficient $b_j$ can be obtained (in principle) by a two-step procedure. First, regress $X_j$ on all the other predictors and take the residuals — this is the part of $X_j$ that is not explained by (uncorrelated with) the rest. Second, do the same for $Y$ — regress $Y$ on the other predictors and take those residuals. Regressing the $Y$-residuals on the $X_j$-residuals gives exactly $b_j$. This is the Frisch–Waugh–Lovell theorem, and it makes precise what "controlling for other variables" actually means: $b_j$ reflects only the covariation between $Y$ and $X_j$ that remains after both have had their shared association with the other predictors stripped out.

Step 1: Specify the model. Choose an outcome $Y$ and a set of predictors believed to influence it, including any suspected confounders. For the coffee example: $Y$ = a cardiovascular risk score, $X_1$ = cups of coffee per day, $X_2$ = cigarettes per day, $X_3$ = hours of exercise per week.

Step 2: Assume linearity and additivity. The model assumes each predictor contributes additively and linearly to $Y$ (this can be relaxed with interaction terms or transformations, but the baseline model is additive-linear). Errors $\varepsilon_i = Y_i - \hat Y_i$ are assumed independent, normally distributed, with constant variance (homoscedasticity) — the same assumptions as simple regression, now applied to a hyperplane.

Step 3: Solve for coefficients via least squares. In matrix form, with $\mathbf{X}$ the design matrix (a column of 1's plus one column per predictor) and $\mathbf{Y}$ the response vector, the least-squares solution is

No new principle is introduced — it is the same minimization as simple regression, generalized to matrix algebra. The $(\mathbf{X}^\top\mathbf{X})^{-1}$ term is where multicollinearity does its damage: when predictors are highly correlated, this matrix becomes nearly singular (non-invertible), and the coefficient estimates become numerically unstable with enormous standard errors.

Step 4: Interpret each $b_j$ as an adjusted effect. The value of $b_1$ (coffee) now represents the change in cardiovascular risk per additional cup of coffee, among people with the same smoking level and same exercise level. If $b_1$ shrinks toward zero once smoking and exercise are added to the model, this is direct numerical evidence that the raw coffee–disease association was substantially confounded.

Step 5: Quantify uncertainty and test hypotheses. Each $b_j$ has a standard error, and $t = b_j / SE(b_j)$ tests $H_0: \beta_j = 0$ (no partial effect of $X_j$, holding others constant) against a $t$-distribution with $n - k - 1$ degrees of freedom. An overall $F$-test assesses whether the full set of predictors explains significantly more variance than a model with none. Nested models can be compared with a partial $F$-test, asking whether adding one or more new predictors significantly improves fit beyond what the simpler model already achieves.

Step 6: Judge overall fit. $R^2$ (coefficient of multiple determination) is the proportion of variance in $Y$ explained by all predictors together. Because $R^2$ mechanically increases whenever any predictor is added — even a useless one — adjusted $R^2$ penalizes for the number of predictors:

This lets researchers compare models with different numbers of predictors fairly.

Multiple regression is a mathematical technique, so its "evidence" is somewhat different from a biological hypothesis: the case for it rests on (1) the formal derivation showing it recovers unbiased partial effects under a set of stated assumptions, and (2) empirical case studies where regression-adjusted findings were later confirmed or refuted by randomized experiments — showing when adjustment succeeds and when it fails.

Framingham Heart Study (begun 1948, Framingham, Massachusetts). This is the canonical biological demonstration of multiple regression's power. Researchers enrolled over 5,000 residents and followed them for decades, measuring blood pressure, cholesterol, smoking, weight, and eventual cardiovascular events. A simple regression of heart disease on cholesterol alone was confounded by age, blood pressure, and smoking — older, hypertensive smokers tend to have different cholesterol profiles too. By fitting a multiple regression (and later, multiple logistic regression, since the outcome was binary — foreshadowing Logistic Regression for Binary Outcomes) with cholesterol, systolic blood pressure, smoking status, and age simultaneously, researchers extracted partial coefficients showing that each of these factors carried an independent, statistically distinguishable risk contribution. This work directly established the concept of "risk factors" as we use it today and demonstrated that adjustment could separate entangled causes in genuinely observational human data — a result later corroborated by randomized antihypertensive and lipid-lowering drug trials that confirmed the directions (though not always the exact magnitudes) of the Framingham partial effects.

Fisher's agricultural field experiments (R.A. Fisher, Rothamsted Experimental Station, 1920s). Although Fisher is most associated with randomization and analysis of variance, his work at Rothamsted also demonstrated the necessity of statistical control when full randomization of every variable is impossible. Crop yield varies with rainfall, soil nitrogen, and sunlight, none of which experimenters can hold fixed across a growing season. Fisher showed that including these measured environmental covariates as additional regressors alongside the treatment variable (fertilizer type) reduced residual variance and sharpened the estimated treatment effect — the seed of the modern method called "analysis of covariance," a direct multiple-regression technique that partitions variance between designed treatments and background confounders measured, not randomized.

The retrospective test case: Hormone Replacement Therapy (HRT) and heart disease. This episode is often cited as a cautionary tale about the limits of adjustment. Large observational cohort studies through the 1980s–90s (e.g., the Nurses' Health Study) used multiple regression to adjust for age, smoking, and socioeconomic status, and consistently found that women on HRT had lower rates of coronary heart disease — suggesting a protective effect. But the Women's Health Initiative (WHI), a large randomized controlled trial (results published 2002), found the opposite: HRT modestly increased cardiovascular risk. The discrepancy arose because regression can only adjust for confounders that are measured and included in the model. HRT users in the observational cohorts were healthier, wealthier, and more health-conscious in ways not fully captured by the covariates used ("healthy user bias") — an unmeasured confounder that no amount of regression could fix. This case is now taught precisely because it shows the boundary of the method: multiple regression adjusts for what you measure, not for what you don't. Together, these episodes justify the technique's use while marking its limits: regression adjustment recovers real, useful partial effects when the relevant confounders are measured (Framingham, Rothamsted), but can be badly misleading when key confounders are omitted (HRT) — a problem no amount of algebra can repair, only design (randomization) can.

Example 1 — A physiological/clinical scenario: birth weight, gestational age, and maternal smoking. Suppose a study of $n = 200$ newborns fits:

where GestAge ranges roughly 32–42 weeks and Smoking is an indicator (1 = mother smoked during pregnancy, 0 = did not). Interpreting each coefficient with the "holding constant" rule:

• $b_0 = 500$: the (extrapolated, not directly meaningful) intercept — predicted weight at 0 weeks gestation and no smoking; not biologically interpretable on its own, just an anchor.
• $b_1 = 180$: for two infants with the same smoking status, each additional week of gestation predicts 180 g more birth weight.
• $b_2 = -210$: for two infants at the same gestational age, a mother who smoked is predicted to have a baby 210 g lighter.

Now suppose a simple regression of birth weight on smoking alone (ignoring gestational age) gives a much larger gap, say $-340$ g. Why the difference? Smokers in the sample deliver slightly earlier on average (say 38.2 weeks vs. 39.5 weeks for non-smokers), and earlier gestational age independently reduces birth weight. The raw $-340$ g estimate is a mix of the true direct effect of smoking on fetal growth and an indirect pathway through shortened gestation. Multiple regression's $-210$ g isolates the effect of smoking that operates independent of gestational age — a genuinely different (and more mechanistically interpretable) quantity. Suppose we compute a specific prediction: a baby born at 36 weeks to a smoking mother:

The exercise demonstrates the mechanics even though real coefficients would differ; the key teaching point is the arithmetic decomposition of a raw vs. adjusted effect via a shared intermediate variable (gestational age).

Example 2 — An ecological/asymmetric-design scenario: fish growth, water temperature, and food density. A biologist measures growth rate ($Y$, mm/month) in a stream fish across 45 sites, along with water temperature ($X_1$, °C) and invertebrate prey density ($X_2$, individuals/m²). The fitted model:

with standard errors $SE(b_1) = 0.22$ and $SE(b_2) = 0.005$, $n = 45$. Compute $t$-statistics for each partial coefficient ($df = n - k - 1 = 45 - 2 - 1 = 42$):

Both exceed the critical value $t_{0.05(2),42} \approx 2.02$, so both temperature and prey density have statistically significant partial effects on growth, each holding the other constant. Now the asymmetry that makes this a genuine test of understanding: suppose temperature and prey density are themselves positively correlated across sites ($r = 0.58$) — warmer sites simply have more invertebrates. A simple regression of growth on temperature alone yields $b_1' = 1.30$ mm/month per °C — noticeably larger than the partial coefficient of 0.85. The difference (1.30 − 0.85 = 0.45) is the portion of temperature's raw association that is actually attributable to its correlation with prey density (warmer sites have more food, and food itself drives growth). The partial coefficient of 0.85 is the direct thermal effect on growth (e.g., via metabolic rate), stripped of the indirect food-mediated pathway. This is exactly the situation multiple regression is built to disentangle — and it is also a caution: if an unmeasured third variable (say, dissolved oxygen) drives both temperature and growth independently, even this adjusted 0.85 could still be biased. Adjustment is only as good as the confounders you actually measured.

Multiple regression does not trace to one single discoverer the way, say, the polymerase chain reaction traces to Kary Mullis. It emerged from the convergence of least-squares theory (Legendre, Gauss, early 1800s), correlation and regression concepts (Francis Galton, 1880s), and the rigorous statistical framework — degrees of freedom, F-tests, experimental design, and the analysis of covariance — built by Ronald Aylmer Fisher at Rothamsted Experimental Station in the 1920s. Fisher's contribution was to turn regression from a descriptive curve-fitting tool into an inferential one: he supplied the sampling theory (t- and F-distributions applied to regression coefficients) that lets a biologist say not just "here is the best-fit hyperplane" but "here is how confident we can be that this partial effect is real." His work on analysis of covariance — blending designed experimental treatments with measured background covariates in a single linear model — is the direct conceptual ancestor of modern multiple regression as used throughout biology and medicine.

Clinical epidemiology and risk prediction. Multivariable regression (linear and, more often, logistic) underlies cardiovascular risk calculators, cancer prognosis scores, and clinical prediction rules — always adjusting simultaneously for age, sex, comorbidities, and the exposure of interest.

Genome-wide association studies (GWAS). Regression of a trait or disease status on millions of genetic markers must adjust for population stratification (ancestry-related confounding) using principal components as covariates, or spurious genotype–phenotype associations swamp the true signal.

Ecology and conservation biology. Species distribution models regress abundance or occurrence on multiple environmental predictors (temperature, precipitation, land use, elevation) simultaneously, disentangling correlated climate and habitat variables to project range shifts under climate change.

Agriculture and animal breeding. Multiple regression (and its extension, mixed-model "animal models") estimates breeding values for traits like milk yield or growth rate while adjusting for herd, year, and management effects — the direct descendant of Fisher's original analysis-of-covariance work.

Public health policy. Studies of air pollution and mortality adjust for age structure, smoking prevalence, and socioeconomic status across cities or time periods, allowing regulators to estimate the health burden attributable specifically to pollutant exposure.

Toward binary outcomes. Many biological outcomes are not continuous but categorical — alive/dead, infected/not, disease/no disease. The linear multiple regression framework extends naturally, with modification of the link function, to Logistic Regression for Binary Outcomes, preserving the same core idea of partial, adjusted effects while modeling probabilities instead of means.

Ask the tutor: How do I choose which confounders to include in my regression model without accidentally adjusting for a mediator and erasing the real effect I'm trying to measure?