Multiple regression and diagnostics

In 5.1 you found that systolic blood pressure rises with age. But older patients also tend to carry more weight, and extra weight raises blood pressure on its own. So part of that age trend might be body mass travelling under the label of age. Multiple regression lets you hold one variable steady while you read the effect of another. It answers a sharper question: what does age do to blood pressure among people of the same build? This lesson fits a model with two predictors, reads the adjusted effects, and then checks whether the fitted line is worth trusting.

The running example is a hypertension screening at a Klinik Kesihatan in Seremban. For twelve adults you recorded age in years, body mass index (BMI) in kg/m², and systolic blood pressure (SBP) in mmHg. The outcome is SBP, still written . Now there are two explanatory variables: age as and BMI as .

The multiple regression model

With one predictor the model was a line. With two it is a plane: each outcome is a sum of an intercept, a term for each predictor, and an error.

The fitting rule has not changed. Least squares picks the intercept and the two slopes that make the sum of squared residuals as small as possible. Each residual is the gap between an observed SBP and the value the plane predicts. What changes is the meaning of each slope.

Adjusted effects and confounding

Fit the simple model first, then add BMI, and watch the age slope move. On its own, age carries a slope of 0.621 mmHg per year. Add BMI and the age slope falls to 0.231. The two predictors are correlated (older patients in this clinic tend to have higher BMI, correlation about 0.79), so the crude age slope was partly BMI working through it. Once you hold BMI fixed, most of the apparent age effect is gone, and what remains is the part age contributes by itself.

age <- c(35, 40, 44, 48, 52, 56, 58, 62, 65, 68, 70, 73)
bmi <- c(24.0, 22.5, 27.0, 25.0, 29.0, 26.0, 31.0, 28.0, 30.0, 27.5, 33.0, 30.0)
sbp <- c(124, 121, 133, 127, 138, 133, 144, 137, 142, 140, 148, 145)

crude <- lm(sbp ~ age)
# now add bmi and compare the age coefficient

crude <- lm(sbp ~ age)
print(round(coef(crude)["age"], 3))   # 0.621 mmHg/year, crude

m <- lm(sbp ~ age + bmi)
print(round(coef(m), 3))              # Intercept 67.666, age 0.231, bmi 1.996
print(round(confint(m), 3))           # 95% CI for each coefficient
# the age slope fell from 0.621 to 0.231 once BMI was held fixed

import pandas as pd
import statsmodels.formula.api as smf
df = pd.DataFrame({"age": [35, 40, 44, 48, 52, 56, 58, 62, 65, 68, 70, 73],
                   "bmi": [24.0, 22.5, 27.0, 25.0, 29.0, 26.0, 31.0, 28.0, 30.0, 27.5, 33.0, 30.0],
                   "sbp": [124, 121, 133, 127, 138, 133, 144, 137, 142, 140, 148, 145]})

crude = smf.ols("sbp ~ age", data=df).fit()
# now add bmi and compare the age coefficient

crude = smf.ols("sbp ~ age", data=df).fit()
print(round(crude.params["age"], 3))   # 0.621 mmHg/year, crude

m = smf.ols("sbp ~ age + bmi", data=df).fit()
print(m.params.round(3).to_dict())     # Intercept 67.67, age 0.231, bmi 1.996
print(m.conf_int().round(3).values)    # 95% CI for each coefficient
# the age slope fell from 0.621 to 0.231 once BMI was held fixed

Reading the fitted plane

The fitted model is SBP = 67.67 + 0.231 × age + 1.996 × BMI. Read each slope holding the other fixed. Among patients of the same BMI, each extra year of age adds about 0.23 mmHg (95% CI 0.12 to 0.34, p = 0.001). Among patients of the same age, each extra unit of BMI adds about 2.0 mmHg (95% CI 1.55 to 2.45, p < 0.001). Both intervals exclude zero, so each predictor carries an independent association after adjusting for the other.

To predict, substitute both values into the equation. Take a 55-year-old with a BMI of 28.

1. Age term: .
2. BMI term: .
3. Add the intercept: mmHg.

The intercept of 67.67 is the predicted SBP when both age and BMI are zero, which is meaningless here. As in 5.1, an intercept far outside the data is a number the model needs, not a reading. Centre both predictors if you want a readable one.

Goodness of fit: R-squared and adjusted R-squared

R-squared still means the share of the outcome's variation the model explains. This model reaches , so age and BMI together account for about 98 percent of the spread in SBP. There is a catch with several predictors: R-squared can only rise when you add a variable, even a useless one, because the fit has more freedom to chase the points. Judging models by R-squared alone always rewards the bigger model.

Adjusted R-squared fixes that by charging a penalty for each predictor. It rises only when a new variable explains more than chance would predict, and it falls when a variable earns its keep.

Here is the sample size and is the number of predictors. Add a column of junk to this model and you can see the two diverge: R-squared edges up from 0.9827 to 0.9836, while adjusted R-squared drops from 0.9789 to 0.9775. The plain figure rewards the noise; the adjusted figure punishes it. Compare models on the adjusted version.

age <- c(35, 40, 44, 48, 52, 56, 58, 62, 65, 68, 70, 73)
bmi <- c(24.0, 22.5, 27.0, 25.0, 29.0, 26.0, 31.0, 28.0, 30.0, 27.5, 33.0, 30.0)
sbp <- c(124, 121, 133, 127, 138, 133, 144, 137, 142, 140, 148, 145)
junk <- c(1, -2, 0, 3, -1, 2, 0, -3, 1, 2, -1, 0)

m <- lm(sbp ~ age + bmi)
# print R-squared and adjusted R-squared, then add junk

m <- lm(sbp ~ age + bmi)
print(round(summary(m)$r.squared, 4))      # 0.9827
print(round(summary(m)$adj.r.squared, 4))  # 0.9789
print(summary(m)$fstatistic[1])            # F = 256 on (2, 9) df

mj <- lm(sbp ~ age + bmi + junk)
print(round(summary(mj)$r.squared, 4))     # 0.9836 , rose
print(round(summary(mj)$adj.r.squared, 4)) # 0.9775 , fell

import pandas as pd
import statsmodels.formula.api as smf
df = pd.DataFrame({"age": [35, 40, 44, 48, 52, 56, 58, 62, 65, 68, 70, 73],
                   "bmi": [24.0, 22.5, 27.0, 25.0, 29.0, 26.0, 31.0, 28.0, 30.0, 27.5, 33.0, 30.0],
                   "sbp": [124, 121, 133, 127, 138, 133, 144, 137, 142, 140, 148, 145]})
df["junk"] = [1, -2, 0, 3, -1, 2, 0, -3, 1, 2, -1, 0]

m = smf.ols("sbp ~ age + bmi", data=df).fit()
# print R-squared and adjusted R-squared, then add junk

m = smf.ols("sbp ~ age + bmi", data=df).fit()
print(round(m.rsquared, 4), round(m.rsquared_adj, 4))   # 0.9827 0.9789
print(round(m.fvalue, 1))                               # 256.0 , F on (2, 9) df

mj = smf.ols("sbp ~ age + bmi + junk", data=df).fit()
print(round(mj.rsquared, 4), round(mj.rsquared_adj, 4)) # 0.9836 0.9775
# R-squared rose; adjusted R-squared fell, flagging the junk

Regression diagnostics: is the model trustworthy?

A model can fit well by R-squared and still rest on broken assumptions. Least squares inference assumes three things about the residuals, and each has a base plot that reveals it.

• Linearity and constant variance. Plot residuals against fitted values. If the model is right, the cloud is a flat, even band around zero with no curve and no fanning. A bend means the straight-line form is wrong. A funnel that widens means the variance grows with the fitted value (heteroscedasticity).
• Normality of residuals. A normal quantile plot (QQ plot) puts the sorted residuals against the values a normal distribution would give. Points on the diagonal mean the residuals are roughly normal; a curve or heavy tails mean they are not.

Run both plots for the Seremban model. The residuals scatter evenly around zero with no clear curve or funnel, and the QQ points sit close to the line, so the assumptions hold here.

age <- c(35, 40, 44, 48, 52, 56, 58, 62, 65, 68, 70, 73)
bmi <- c(24.0, 22.5, 27.0, 25.0, 29.0, 26.0, 31.0, 28.0, 30.0, 27.5, 33.0, 30.0)
sbp <- c(124, 121, 133, 127, 138, 133, 144, 137, 142, 140, 148, 145)
m <- lm(sbp ~ age + bmi)

# plot residuals vs fitted, then a normal QQ plot

par(mfrow = c(1, 2))
plot(fitted(m), residuals(m), pch = 19,
     xlab = "Fitted SBP", ylab = "Residual",
     main = "Residuals vs fitted")
abline(h = 0, col = "grey", lty = 2)

qqnorm(residuals(m), pch = 19)   # normal quantile plot
qqline(residuals(m), col = "#B22222")
# even band around zero, points near the line: assumptions hold

import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import statsmodels.formula.api as smf
df = pd.DataFrame({"age": [35, 40, 44, 48, 52, 56, 58, 62, 65, 68, 70, 73],
                   "bmi": [24.0, 22.5, 27.0, 25.0, 29.0, 26.0, 31.0, 28.0, 30.0, 27.5, 33.0, 30.0],
                   "sbp": [124, 121, 133, 127, 138, 133, 144, 137, 142, 140, 148, 145]})
m = smf.ols("sbp ~ age + bmi", data=df).fit()

# plot residuals vs fitted, then a normal QQ plot

fig, ax = plt.subplots(1, 2, figsize=(9, 4))
ax[0].scatter(m.fittedvalues, m.resid, color="black")
ax[0].axhline(0, color="grey", linestyle="--")
ax[0].set_xlabel("Fitted SBP"); ax[0].set_ylabel("Residual")
ax[0].set_title("Residuals vs fitted")

stats.probplot(m.resid, plot=ax[1])   # normal QQ plot
plt.show()

print(round(stats.shapiro(m.resid).pvalue, 3))  # 0.503 , no evidence against normality

Influential points and outliers

A single odd observation can pull the whole plane toward it. Two ideas separate the harmless from the dangerous. Leverage measures how far a patient's predictor values sit from the average, so a patient with an unusual age-and-BMI combination has high leverage. Cook's distance combines leverage with the size of the residual to measure how much the fitted coefficients would move if you dropped that patient. A point can have a large residual yet little influence, or sit near the line yet swing the model. Cook's distance catches the ones that matter.

Suppose a clerk enters a thirteenth patient, a 50-year-old, but types the BMI as 38.0 instead of the true 28.0, with SBP 130. Fit the model and read the influence measures. Cook's distance for that patient is 10.37, against 0.17 or less for everyone else. Its leverage is 0.77, far above the others, because a BMI of 38 sits well outside the rest of the cloud. A common rule of thumb flags Cook's distance above , here , and this point clears that by thirtyfold.

1. Fit the model and pull out Cook's distance and leverage for every patient.
2. Flag any point with Cook's distance above , or leverage well above the average .
3. Inspect the flagged record. Here BMI 38 is a data-entry error, not a real patient.
4. Refit without it (or after fixing the typo) and check whether the conclusions hold.

Refitting without the bad row returns the age slope to 0.231 and the BMI slope to 1.996. With the error left in, those slopes distort to 0.573 and 0.346, which would have told a clinician the wrong story about both predictors.

age <- c(35, 40, 44, 48, 52, 56, 58, 62, 65, 68, 70, 73, 50)
bmi <- c(24.0, 22.5, 27.0, 25.0, 29.0, 26.0, 31.0, 28.0, 30.0, 27.5, 33.0, 30.0, 38.0)
sbp <- c(124, 121, 133, 127, 138, 133, 144, 137, 142, 140, 148, 145, 130)
m <- lm(sbp ~ age + bmi)

ck <- cooks.distance(m)
# find the worst point, then refit without it

ck <- cooks.distance(m)
hv <- hatvalues(m)
print(round(ck, 2))                  # patient 13 = 10.37 , rest below 0.17
print(round(hv, 2))                  # patient 13 = 0.77 , high leverage
print(which(ck > 4 / length(sbp)))   # flags patient 13

m2 <- lm(sbp[-13] ~ age[-13] + bmi[-13])
print(round(coef(m2), 3))            # age 0.231, bmi 1.996 restored

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
df = pd.DataFrame({"age": [35, 40, 44, 48, 52, 56, 58, 62, 65, 68, 70, 73, 50],
                   "bmi": [24.0, 22.5, 27.0, 25.0, 29.0, 26.0, 31.0, 28.0, 30.0, 27.5, 33.0, 30.0, 38.0],
                   "sbp": [124, 121, 133, 127, 138, 133, 144, 137, 142, 140, 148, 145, 130]})
m = smf.ols("sbp ~ age + bmi", data=df).fit()

infl = m.get_influence()
# read Cook's distance and leverage, then refit without the worst point

infl = m.get_influence()
ck = infl.cooks_distance[0]
hat = infl.hat_matrix_diag
print(np.round(ck, 2))               # patient 13 = 10.37 , rest below 0.17
print(np.round(hat, 2))              # patient 13 = 0.77 , high leverage
print(np.where(ck > 4 / len(df))[0] + 1)  # flags patient 13

clean = df.drop(index=12)
m2 = smf.ols("sbp ~ age + bmi", data=clean).fit()
print(m2.params.round(3).to_dict())  # age 0.231, bmi 1.996 restored

Common mistakes

• Reading a partial slope as a marginal one. The adjusted age slope of 0.231 is the effect among people of the same BMI, not the total age effect you would see across the clinic. Drop the "holding BMI fixed" and you are quoting the wrong number.
• Comparing models by R-squared. R-squared never falls when you add a predictor, so it always favours the larger model. Use adjusted R-squared, or a measure like AIC, when the models differ in size.
• Skipping the residual plots. A high R-squared says nothing about whether the form is right. A curved residuals-versus-fitted plot means a straight model on a bent relationship, and no fit statistic warns you.
• Deleting outliers on sight. A point with high influence is a prompt to check the record, not a licence to drop it. Remove it only when it is a genuine error, and report results with and without it.
• Confusing a large residual with high influence. A far-off point near the centre of the predictors may barely move the fit, while a modest residual at an extreme predictor value can swing it. Cook's distance, not the residual alone, tells you which.

Tips

• Show the crude and adjusted slope side by side. The reader learns more from "age slope 0.62 crude, 0.23 adjusted for BMI" than from either number alone. The shift is the confounding story.
• Name what you adjusted for, every time. An adjusted effect is only adjusted for the variables in the model. Write "adjusted for age and BMI" so no one reads it as fully controlled.
• Make the diagnostic plots part of the fit, not an afterthought. Residuals-versus-fitted and a QQ plot take two lines and catch most problems before they reach a report.
• Treat Cook's distance as a triage tool. Sort by it, inspect the top few records against the source data, and act only on real errors.
• Watch for collinearity when both slopes look unstable. If two predictors are strongly correlated and both confidence intervals are wide, the data cannot separate them cleanly. Consider dropping one or combining them.