How do several factors together predict an outcome?

Start here

What does "holding the other predictors constant" mean?

Multiple linear regression gives each predictor its own coefficient. Each coefficient is the difference in the outcome for a 1-unit difference in that predictor, among patients who have the same values of every other predictor in the model.

The right-hand panel shows how the model works. The coefficient for age is estimated only from the variation in age that the other predictors do not explain. Two patients of different ages who have had diabetes for the same time are what inform the age coefficient.

y = b₀ + b₁x₁ + b₂x₂ + … + bₖxₖ

HbA1c = 7.208 + 0.004 × age + 0.080 × duration

y

the outcome (dependent variable), here HbA1c in %

x₁ … xₖ

the k predictors (independent variables)

b₀

the intercept: predicted outcome when every predictor is 0, often a value outside the data (age 0)

b₁ … bₖ

the partial regression coefficients: each one is adjusted for the other predictors in the model

In this model, age has B = 0.004 (95% CI −0.019 to 0.027, p = 0.73) and duration has B = 0.080 (0.035 to 0.125, p < 0.001). The coefficient is adjusted only for the predictors in the model, as they were measured. A confounder that was not measured is not adjusted for.

Interactive

Add or remove predictors and watch the coefficients, R² and adjusted R² change

Goodness of fit

How much of the variation does the model explain?

The total variation in the outcome is split into the part the model explains and the part left over. R² is the explained share.

SST = Σ(y − ȳ)²   SSR = Σ(ŷ − ȳ)²   SSE = Σ(y − ŷ)²   SST = SSR + SSE

R² = SSR ÷ SST = 1 − SSE ÷ SST    R = √R²

adjusted R² = 1 − (1 − R²) × (n − 1) ÷ (n − k − 1)

ŷ

the value the model predicts for each patient

SSR

regression (explained) sum of squares

SSE

residual (error) sum of squares: what the model leaves unexplained

R

multiple correlation coefficient: the correlation between the observed and the predicted outcome. With one predictor it equals the size of Pearson r

n, k

number of patients, number of predictor terms (not counting the intercept)

Why R² always rises. Adding a predictor can never make the least squares fit worse, because the new coefficient could be set to 0 and the fit would be the same as before. So R² goes up, or stays level, even for a column of random numbers. Adjusted R² goes up only when the added predictor's t statistic is larger than 1 in absolute value, and it can be negative.

Multiple R for the full teaching model is √0.513 = 0.716: the observed HbA1c values correlate at 0.716 with the values the model predicts.

Interactive

Add useless predictors and watch R² rise while adjusted R² does not

Tests

What do the F test and the t tests each test?

The F test asks whether the predictors together explain more than chance would. Each t test asks whether one coefficient differs from 0 when the other predictors are held constant.

F = (SSR ÷ k) ÷ [SSE ÷ (n − k − 1)], with k and n − k − 1 degrees of freedom

t = B ÷ SE(B), with n − k − 1 degrees of freedom    95% CI = B ± t₀.₉₇₅ × SE(B)

• Overall F: tests whether all slope coefficients are 0 at once. Here F(5, 114) = 24.05, p < 0.001. A significant F says the model explains something, not that every predictor matters.
• t for each coefficient: duration t = 5.11, p < 0.001. Age t = 0.01, p = 0.99: once duration, BMI and activity are known, age adds nothing.
• F for a block of terms: compares two nested models, the full model and the model without the block. Physical activity (two 0/1 variables) gave F(2, 114) = 12.99, p < 0.001, and added 0.111 to R².
• The intercept's p value tests whether the predicted outcome is 0 when every predictor is 0. For HbA1c at age 0 this has no clinical meaning.

A significant overall F with no significant t test is a sign that two or more predictors overlap. That is covered in the multicollinearity chapter.

Partial correlation

How strongly is one predictor related to the outcome once the others are allowed for?

Partial correlation removes the effect of the other predictors from both the outcome and the predictor. Semi-partial (part) correlation removes it from the predictor only.

partial r = t ÷ √(t² + df)    semi-partial r² = t² × (1 − R²) ÷ df

For duration in the full model, t = 5.11 and df = 114: partial r = 0.43 and semi-partial r = 0.33 (semi-partial r² = 0.111). For BMI, partial r = 0.39 and semi-partial r = 0.29. The partial correlation is always at least as large as the semi-partial one, because it divides by a smaller amount of variation.

Standardised coefficients

Can you compare coefficients measured in different units?

B is in outcome units per unit of the predictor. The standardised coefficient β is in standard deviations of the outcome per standard deviation of the predictor.

β = B × SD(x) ÷ SD(y)

Duration: β = 0.0959 × 6.075 ÷ 1.046 = 0.56. BMI: β = 0.31. Sedentary: β = 0.37. A patient whose duration is 1 SD (6.1 years) longer has HbA1c about 0.56 SD (0.58 percentage points) higher, holding the other predictors constant.

• β depends on the SDs in your sample. The same biology gives a different β in a sample with a narrower range of duration.
• For a 0/1 predictor, a change of 1 SD has no clinical meaning. Report B for categories.
• Report B with its units and 95% CI as the main result. β can be added, but a table of β alone cannot be read in clinical units.

Categorical predictors

How do you put a predictor with several categories into the model?

A predictor with k categories becomes k − 1 indicator (dummy) variables, each coded 0 or 1. The category with no variable of its own is the reference category.

• Why k − 1: a variable for every category would add up to 1 for every patient, the same as the intercept column, and the coefficients could not be estimated.
• Choosing the reference: a clinically natural comparison group with enough patients. A reference category with very few patients gives wide intervals for every comparison.
• Testing the predictor: use one F test for the whole set of dummies. A single non-significant dummy does not mean that category can be merged or dropped.
• Never enter codes 1, 2, 3 as one number for unordered categories. That forces the same step between each pair of categories.

Interaction terms

What if the effect of one predictor depends on another?

An interaction term is the product of two predictors. Its coefficient is how much the slope of one predictor changes for a 1-unit change in the other.

HbA1c = b₀ + … + b₁ × (BMI − 27) + b₂ × sedentary + b₃ × (BMI − 27) × sedentary

• Main effects change meaning. With the interaction in the model, the sedentary coefficient (0.85) is the difference at BMI 27, and the BMI coefficient is the slope in patients who are not sedentary.
• Centre the continuous predictor (subtract a sensible value such as 27) so that the main effects describe a realistic patient. Without centring, the sedentary coefficient is the difference at BMI 0: −2.16 in the same data.
• Test the interaction term itself. A significant effect in one group and not in the other does not show that the effects differ.
• Plan interactions before the analysis. Testing every possible product term produces chance findings.

Assumptions

How do you check that the model is suitable for the data?

The assumptions are about the residuals, the differences between the observed and the predicted outcome. Plot them before you read the p values.

1. Linearity: the outcome changes in a straight line with each continuous predictor, holding the others constant. Check: residuals against fitted values and against each predictor show no curve.
2. Independence: one row per patient, with no repeated measurements or clusters analysed as if independent.
3. Constant variance: the spread of the residuals is similar at all fitted values. Check: no funnel shape.
4. Normal residuals: needed for accurate p values and confidence intervals, most of all in small samples. Check: a histogram or normal Q-Q plot of the residuals. The raw outcome does not need to be normal.
5. No single patient controls the result: check outliers and influential points (for example Cook's distance), and the case records behind them.
6. Predictors that do not overlap too much: check variance inflation factors (multicollinearity chapter).

Interactive

Change how the data were generated and watch the residual plot

Decide

How should each predictor enter the model?

Answer the questions for each predictor before you run the analysis.

Comparison

Three pairs that residents confuse

Reading a paper

How do you read a regression table in a paper?

Read the table from the rows to the footer: what each row compares, its size in units, its precision, and the model it came from.

n = 120 · R² = 0.513 · adjusted R² = 0.492 · F(5, 114) = 24.05, p < 0.001 · outcome HbA1c (%) · teaching data

1. Outcome and units: find them in the title or footer. Here every B is in percentage points of HbA1c.
2. Each row: a continuous predictor is per unit (per year, per kg/m²). A category row names its reference (“vs active”). If the reference is not stated, the row cannot be read.
3. Size before significance: sedentary patients have HbA1c 0.93 percentage points higher than active patients with the same age, duration and BMI.
4. Precision: the 95% CI for sedentary runs from 0.57 to 1.29. The CI for age includes 0 and is narrow, so a clinically important age effect at the same duration is unlikely in these data.
5. Adjustment set: every coefficient is adjusted for the other rows, and only those. Check whether an obvious confounder is missing, and whether a variable on the causal path was included.
6. The footer: n, R², adjusted R² and the F test. Divide n by the number of predictor terms (here 120 ÷ 5 = 24).
7. What is missing: unadjusted estimates, how predictors were chosen, residual checks and variance inflation factors.

The formula

Worked example: the full HbA1c model

Teaching data, made up: 120 adults with type 2 diabetes. HbA1c on age, duration of diabetes, BMI and physical activity (reference: active). Every number was computed with statsmodels OLS in Python.

1. Sums of squares from the software: SST = 130.15, SSR = 66.81, SSE = 63.33. n = 120, k = 5, residual df = 120 − 5 − 1 = 114.
2. R² = 66.81 ÷ 130.15 = 0.513. Multiple R = √0.513 = 0.716.
3. Adjusted R² = 1 − (1 − 0.513) × 119 ÷ 114 = 0.492.
4. F = (66.81 ÷ 5) ÷ (63.33 ÷ 114) = 13.36 ÷ 0.556 = 24.05, with 5 and 114 df, p < 0.001.
5. Duration: B = 0.0959, SE = 0.0188, t = 5.11 (from the unrounded values), df 114, p < 0.001.
6. 95% CI = 0.0959 ± 1.981 × 0.0188 = 0.059 to 0.133 percentage points per year.
7. Standardised β = 0.0959 × 6.075 ÷ 1.046 = 0.56 (SD of duration 6.075 years, SD of HbA1c 1.046%).
8. Partial r = 5.11 ÷ √(5.11² + 114) = 0.43. Semi-partial r² = 5.11² × (1 − 0.513) ÷ 114 = 0.111, the fall in R² if duration is removed.
9. Prediction for a 60-year-old with diabetes for 10 years, BMI 30, sedentary: 4.777 + 0.0001 × 60 + 0.0959 × 10 + 0.0836 × 30 + 0.932 = 9.18%. 95% CI for the mean of such patients 8.89 to 9.48. 95% prediction interval for one such patient 7.68 to 10.69.

import statsmodels.formula.api as smf
# d has one row per patient. moderate and sedentary are 0/1 columns (reference: active)
fit = smf.ols("hba1c ~ age + duration + bmi + moderate + sedentary", data=d).fit()
print(fit.summary())                      # B, SE, t, p, 95% CI, R-squared, F
print(fit.rsquared, fit.rsquared_adj)     # 0.513, 0.492

In a real paper

Multiple linear regression in published studies

How to run it

Running multiple linear regression in SPSS, R and jamovi

Each package fits the same model. Enter categorical predictors as indicator variables or as factors, and ask for confidence intervals and residual plots.

IBM SPSS Statistics

1. Analyze > Regression > Linear... Select the numeric dependent variable and one or more numeric independent variables. Independent variables can be grouped into blocks, each with its own entry method.
2. Leave Method as Enter, which enters all variables in a block in a single step.
3. The independent variables must be numeric, so enter physical activity as two 0/1 variables you have created, leaving out the reference category.
4. Statistics: Estimates gives B, its standard error, the standardized coefficient beta, t and the two-tailed significance. Also tick Confidence intervals, Model fit (multiple R, R², adjusted R², standard error of the estimate and the ANOVA table), R squared change for blocks, and the part and partial correlations.
5. Plots: put *ZRESID on the Y axis and *ZPRED on the X axis to check linearity and equal variance. Tick the histogram and the normal probability plot of the standardized residuals. “Produce all partial plots” shows each predictor after the others are removed.
6. Partial correlation on its own: Analyze > Correlate > Partial..., with the other predictors as control variables.

R

d$activity <- relevel(factor(d$activity), ref = "active")   # choose the reference category

fit <- lm(hba1c ~ age + duration + bmi + activity, data = d)  # treatment contrasts: k - 1 dummies against the reference
summary(fit)      # B, SE, t, p, R-squared, adjusted R-squared, F statistic
confint(fit)      # 95% confidence intervals

fit_no_act <- lm(hba1c ~ age + duration + bmi, data = d)
anova(fit_no_act, fit)   # F test for the activity block (nested models)

d$bmi_c <- d$bmi - 27
fit_int <- lm(hba1c ~ age + duration + moderate + bmi_c * sedentary, data = d)  # 0/1 columns, and a*b = a + b + a:b

plot(fit)                                   # residuals vs fitted, Q-Q, scale-location, residuals vs leverage
lm.beta::lm.beta(fit)                       # standardised coefficients (package lm.beta)
ppcor::pcor(d[, c("hba1c", "age", "duration", "bmi")])    # partial correlations (package ppcor)
ppcor::spcor(d[, c("hba1c", "age", "duration", "bmi")])   # semi-partial correlations

jamovi

1. Analyses → Regression → Linear Regression.
2. Move HbA1c to Dependent Variable, continuous predictors to Covariates, and categorical predictors to Factors. The reference level of each factor can be set (the refLevels option of linReg in the jmv R package).
3. Model Fit: tick R, R² and Adjusted R². Model Coefficients: tick ANOVA test, Confidence interval and Standardized estimate.
4. Assumption Checks: tick Q-Q plot of residuals, Residuals plots and Cook's distance.
5. To compare nested models, use Model Builder and add a second block with “+ Add New Block”.
6. Partial and semi-partial correlations are in the Partial Correlation analysis (corrPart in the jmv R package, with type partial or semipartial).

How to report it

What to write in the methods and results

Report the model so that a reader could fit it again: the outcome, every predictor and how it was coded, how predictors were chosen, the checks, and each coefficient with its units and 95% CI.

• Methods: the outcome, each predictor with its units or categories and reference category, and why each was included. Any interaction or transformation planned. How assumptions and collinearity were checked, how missing data were handled, and the software and version.
• Results: n in the model, and for each predictor B, 95% CI and p, preferably in a table. The unadjusted estimate for the main exposure. R² or adjusted R², and any interaction tested, with its estimate.

If you enter the residual checks, state what they showed, not only that they were done. The guidance these items come from is in the Guidelines section below.

By specialty

Multiple regression problems by specialty

Common myths

Common mistaken beliefs

• ✗Myth: The coefficient for age is the effect of age.

  ✓It is the difference in the outcome per year of age among patients with the same values of the other predictors in the model. In the teaching data it was 0.036 alone and 0.000 with duration, BMI and activity in the model.

• ✗Myth: A higher R² means a better model.

  ✓R² never falls when a predictor is added, even a column of random numbers. Compare models with adjusted R², and judge them by their purpose and their residual plots.

• ✗Myth: If one category of a factor is not significant, merge or drop it.

  ✓Each dummy coefficient compares one category with the reference. Test the factor as a whole: activity gave F(2, 114) = 12.99, p < 0.001, although moderate activity alone had p = 0.13.

• ✗Myth: The predictor with the largest coefficient matters most.

  ✓B depends on units. HbA1c changes 0.096 per year of diabetes and 0.084 per kg/m² of BMI, and these cannot be ranked without allowing for how much each predictor varies.

• ✗Myth: "Adjusted for" means confounding has been removed.

  ✓The model adjusts only for the variables it contains, as they were measured. Confounders that were not measured, or were measured poorly, can still bias the coefficient.

• ✗Myth: The outcome must be normally distributed.

  ✓The assumption concerns the residuals, the differences between observed and predicted values, not the raw outcome.

Why this matters to you

Five consequences of misreading multiple regression

1. A conclusion you cannot defend at the viva

   Writing that older patients have poorer control when, at the same duration of diabetes, age made no difference.

2. A protocol sent back

   "Multivariate analysis will be done" without naming the outcome, the predictors, their coding and why each was chosen.

3. Wrong numbers from a coding error

   A three-category predictor entered as 1, 2, 3 forces equal steps between categories and gives a coefficient that describes no real comparison.

4. Misreading papers

   Taking a standardised β of −0.04 as if it were in clinical units, or reading a dummy coefficient without knowing its reference category.

5. Predictions that fail for patients

   Using the equation for patients outside the range of the data, or trusting a high R² from a model whose residuals curve.

History

History of multiple regression

1. 1886

   Francis Galton

   Published “Regression towards mediocrity in hereditary stature”: children of unusually tall or short parents were, on average, closer to the population average.

2. 1896

   Karl Pearson

   Published a mathematical treatment of correlation and regression in the Philosophical Transactions of the Royal Society, calculating both from the product-moment.

3. 1922

   R. A. Fisher

   Published “The goodness of fit of regression formulae, and the distribution of regression coefficients” in the Journal of the Royal Statistical Society.

4. 1977

   R. Dennis Cook

   Published “Detection of influential observation in linear regression” in Technometrics.

5. 2013

   Bertha Hidalgo and Michael Goodman

   In one year of the American Journal of Public Health, 25 of 30 articles that described a "multivariate" method had used multivariable analyses with a single outcome.

6. 2015

   Peter Austin and Ewout Steyerberg

   Simulations showed that linear regression needs only about two subjects per variable for accurate coefficients, standard errors and confidence intervals, but many more for an unbiased R². Adjusted R² behaved well.

7. 2019

   Richard Riley and colleagues

   Proposed four criteria for the minimum sample size of a linear regression prediction model, including a global shrinkage factor of 0.9 or more and a difference of 0.05 or less between apparent and adjusted R².

8. 2024

   Gary Collins and colleagues

   Published TRIPOD+AI, a 27-item checklist for reporting prediction models built with regression or machine learning, replacing TRIPOD 2015.

Guidelines

Multiple regression in reporting guidelines

Reading a paper

Five checks on any multiple regression in a paper

• Are the outcome, every predictor, its units and the reference category of each categorical predictor stated?
• Are coefficients given in clinical units with 95% confidence intervals, not only standardised β or p values?
• Were the predictors chosen before the analysis for stated reasons, or by a stepwise procedure?
• Are unadjusted estimates shown alongside adjusted ones, and is the number of patients large for the number of predictor terms?
• Did the authors check residuals, interactions and collinearity, and report what they found?

Why is it called that?

Word surgery

Exam answers

Short-note outlines

Further reading

Related Stateazy articles