AdvancedIndependently reviewed, not yet spot-checked by a human

Regression model diagnostics

What each of the four residual plots is looking for, why passing all four still leaves problems undetected, where the importance of VIF is overstated, how influence is measured, and why stepwise regression is one of the most common statistical errors in the medical literature — demonstrated by simulation.

What diagnostics are checking

Getting a set of coefficients out does not mean the model is right. Every p-value and every confidence interval in linear regression rests on four assumptions:

Correct specification — the shape of the relationship between the variables you put in and the outcome is right (linearity, interactions)
Constant variance (homoscedasticity) — the spread of the residuals does not change with the fitted value
Approximately normal residuals — this mainly affects intervals in small samples; in large samples the central limit theorem rescues you
Independent observations — this one is invisible in the data and has to come from the design (repeated measures, cluster sampling, and patients from the same hospital all violate it). For what to do once it is violated, see repeated measures and clustered data

Diagnostics work through the first three, plus two problems the output never shows: collinearity and influential points. This page reuses the linear model from B2-01 and the logistic model from B2-02, whose output you have already read.

What to look for in the residual plots

Four diagnostic plots: residuals against fitted values, a normal Q-Q plot, a scale-location plot, and leverage against standardised residuals with Cook's distance contours. None of the four shows an obvious systematic pattern. — The four standard diagnostic plots for the B2-01 model. The red line is a lowess smooth; the dashed red curves in panel D are the Cook's distance contours at 0.5 and 1, and the circled point is this model's largest Cook's D.Plotting script `figures/scripts/B2-04-model-diagnostics.R`

Plot	What it looks for	What broken looks like	This model
A. Residuals vs fitted	Whether the specification is right (missing curvature, interactions)	The smooth bends into a U or an inverted U	No overall curvature detected (Tukey test p 0.632, not statistically significant)
B. Normal Q-Q	The tails of the residual distribution	Both ends leave the diagonal (heavy tails, skew)	No departure from normality detected (Shapiro-Wilk p 0.765, not statistically significant)
C. Scale-location	Whether the variance is constant	The smooth climbs to the right (a fan shape)	No change in variance with fitted value detected (Breusch-Pagan p 0.765, not statistically significant)
D. Leverage vs residual	Whether a single point is driving the model	A point falls outside the Cook’s D contour at 0.5 or 1	Largest Cook’s D 0.118

The residual standard deviation over the lower half of the fitted values is 628 and over the upper half 643 — almost identical, so panel C being flat has numbers behind it.

Passing all four plots does not mean the specification is correct.

B2-01 already found that age is not linear in this dataset: adding a quadratic term for age improves the model detectably (p 0.018). Yet the overall curvature test in column A above returns p 0.632 and catches nothing.

The reason is that the horizontal axis of panel A is the fitted value, and the fitted value is a weighted blend of every variable. Curvature in a single variable is diluted inside that blend. Plot the residuals against each continuous variable separately and it shows up:

Residuals plotted against	p for the quadratic term	Caught?
Fitted values (Tukey)	0.632	Missed it
Age	0.017	Caught it
Maternal weight	0.749	No departure detected

So the four standard plots are a starting point, not the finish. Every continuous variable in the model deserves its own residual plot — or use restricted cubic splines and let the data choose the shape.

Two residual plots: on the left, residuals against maternal age, where the lowess smooth turns upward at both ends; on the right, residuals against maternal weight, where the smooth is roughly flat. — The same residuals, plotted against one variable at a time. The smooth in the age panel turns up at both ends — exactly what the overall test in the previous section missed. The maternal weight panel is the control, and looks the way it should.Plotting script `figures/scripts/B2-04-model-diagnostics.R`

Collinearity and VIF

Collinearity means one variable can be very nearly predicted from a linear combination of the others. The variance inflation factor (VIF) measures precisely that:

\mathrm{VIF}_j = \frac{1}{1 - R_j^2}

where $R_j^2$ is the R² from regressing the $j$ -th variable on all the others. $\mathrm{VIF} = 4$ means that coefficient’s standard error has been inflated two-fold.

Every VIF in the B2-01 model is low (the largest is 1.34). To demonstrate the problem, a column holding the same maternal weight converted to kilograms is added deliberately — something that happens often once real datasets are merged. Note that the conversion carries a little measurement / rounding error (lwt * 0.4536 + rnorm(n, 0, 0.5)), which is what merged real data actually looks like; had the conversion been exactly proportional, R would declare the column aliased with lwt, return NA for the coefficient, report an infinite VIF, and the table below could not be produced at all:

Variable	VIF, original model	VIF after adding the duplicate column
age	1.10	1.11
lwt	1.22	771.05
lwt_kg	—	771.66
smoke_fYes	1.16	1.16
race_fBlack	1.19	1.19
race_fOther	1.34	1.35
ht	1.08	1.08
ui	1.04	1.04

Collinearity damages only the variables caught up in it; the other coefficients are untouched. This is where the importance of VIF is most often overstated.

Look at the other numbers from the same refit:

Variable	Coefficient, original	After the duplicate column	SE, original	After
smoke_fYes	-360.7	-358.5	104.0	104.1
ht	-590.0	-585.3	200.3	200.4
ui	-528.5	-520.4	135.1	135.4

The coefficients and standard errors for smoking, hypertension and uterine irritability barely move; the residual standard deviation goes from 647.3 to 647.4, and R² from 0.241 to 0.245.

The only thing wrecked is the variable that was caught up in it: the coefficient for maternal weight flips from 4.40 (p 0.011) to -37.22 (p 0.387), with the standard error inflated 25-fold and even the direction reversed.

What this means in practice:

If the variable you care about is not part of the collinearity, a high VIF is not a problem and there is no need to drop anything for it. When a reviewer asks you to “deal with the collinearity”, first ask which variables are collinear and whether yours is among them.
Collinearity causes no bias; it only inflates variance. The point estimates remain unbiased, just imprecise.
There is no magic VIF threshold. Five and ten are both conventions. The question worth asking is whether this standard error is still small enough to answer your question.
Only two situations genuinely need action: the same thing has been measured twice (as in the demonstration above — drop one), or the collinear variables are only meaningful together (use a summary score or dimension reduction).
⚠️ The VIFs on this page are computed from the definition (each column of the model matrix regressed on the remaining columns), so every dummy column of a multi-category factor gets its own VIF. car::vif() reports a generalised VIF (GVIF) for factors, one value per factor, and the numbers differ from these — know which one you are reading.

Influential points

Outlier, leverage and influence are three different things:

Outlier: $y$ is far from its fitted value → look at the standardised residual
Leverage: $x$ sits at the edge of the covariate space → look at the hat value; the average is $p/n$ = 0.042
Influence: would the coefficients change if the point were removed → Cook’s distance, which combines the other two multiplicatively: it is large only when the residual and the leverage are both large

D_i = \frac{\sum_j (\hat{y}_j - \hat{y}_{j(i)})^2}{p \cdot \hat{\sigma}^2}

Only the third one really matters. A point with an extreme $x$ that sits on the regression line has high leverage and low influence; it is in fact helping to hold the slope steady.

The five largest Cook’s D values in this model:

Row	Birth weight (g)	Maternal age	Maternal weight (lb)	Leverage	Standardised residual	Cook’s D
130	4990	45	123	0.106	2.82	0.118
133	1135	34	187	0.153	-1.71	0.066
132	1021	29	130	0.057	-2.90	0.063
106	3790	25	241	0.143	1.66	0.058
102	3756	19	184	0.108	1.72	0.045

No influential point was detected in this model, and saying so honestly beats manufacturing one for the demonstration.

Three thresholds are in common use: Cook’s D > 1 (the most permissive), > 0.5, and > 4/n (here 0.021). The largest value in this model is 0.118: 0 points exceed 1, 0 exceed 0.5, and 14 exceed 4/n.

The 4/n line flags a crowd of points in any decent-sized sample — it means “relatively large”, not “a problem”. The way to check is to refit: drop the patient with the largest Cook’s D and the smoking coefficient moves from -360.7 to -341.7 (a shift of 5.3%), while hypertension moves from -590.0 to -593.4. No conclusion turns on that one patient.

When you do find a genuinely influential point, what to do and what not to do:

✅ First check whether it is a data error (wrong units, a typo, someone who should never have been included). If it is, correct or exclude it and say so in the paper.
✅ If it is not an error, run a sensitivity analysis: report the result with and without that person, side by side.
❌ Do not delete it because it makes the result look bad. Removing influential observations without declaring it is data manipulation.
❌ Do not use Cook’s D as an automated cleaning rule.

These diagnostics also close a loose end from B2-01. The quadratic age term on that page was statistically supported (p 0.018), but refitting without the 3 mothers older than 35 makes the evidence disappear (p 0.373). 3 people can decide a conclusion in a model of 189 — which is why “significant” and “robust” are two different things, and why splines beat polynomials.

Diagnostics differ for binary and count outcomes

Those four plots cannot be carried over to logistic regression. The reason is simple: the outcome takes only the values 0 and 1, so the raw residuals can only fall on two curves, and the plot is always two bands with nothing to read in it (in this model the raw residuals range from -0.74 to 0.91, but the shape is determined entirely by the outcome value).

Three other things take their place:

1. Calibration — whether the predicted probabilities match the observed rates. Split the patients into ten groups by predicted probability:

Decile	n	Mean predicted probability	Expected events	Observed events
1	19	0.057	1.1	0
2	19	0.102	1.9	2
3	19	0.154	2.9	5
4	19	0.207	3.9	3
5	19	0.243	4.6	4
6	18	0.283	5.1	7
7	19	0.341	6.5	6
8	19	0.450	8.6	8
9	19	0.551	10.5	11
10	19	0.731	13.9	13

The Hosmer-Lemeshow statistic is 4.68 on df = 8, p 0.792 — no systematic miscalibration was detected.

2. Discrimination — see ROC and AUC (the B4 family). Calibration and discrimination are independent: a model can rank patients well while its probabilities sit systematically too high, and it can produce well-calibrated probabilities while ranking poorly.

3. Influence still matters — a glm object has cooks.distance() and hatvalues() as well, and they are read the same way.

As for Poisson regression, the crucial diagnostic is not a residual plot but the dispersion: whether the Pearson chi-square divided by the degrees of freedom is close to 1. On that page the example comes out at 11.05, meaning the standard errors need to be inflated by roughly 3.32-fold.

Stepwise regression: one of the most common statistical errors in the medical literature

Stepwise regression means letting an algorithm add and drop variables automatically by p-value or AIC until a set of “significant” variables remains. R’s step() / stepAIC() and the forward / backward options in SPSS all do this.

It is extremely common in the literature and almost always wrong. A simulation makes it visible: take the B2-01 model, add 10 pure standard normal random variables (unrelated to birth weight in any way), run stepAIC(), and repeat 300 times.

300 simulation runs. Left: stepAIC retains 1.81 pure noise variables on average. Right: of those, 0.63 on average reach p < 0.05 in the final model, and 48% of runs produce at least one "significant" noise variable.Plotting script figures/scripts/B2-04-model-diagnostics.R

A harsher version: replace the outcome with random numbers too, and feed in 15 random predictors, so there is nothing in the data to find at all. Across 300 runs:

2.59 variables retained on average
55% of final models contain at least one variable with p < 0.05
the mean R² is 0.047
64% of final models also return a “significant” overall F test

That last figure needs care with its denominator: in 24 of those 300 runs stepAIC deleted every variable and left only the intercept, and such a model has no overall F test. The 64% above counts those 24 runs as “not significant” (denominator 300); count only the 276 runs that actually have an F test and the figure is 70%. Both numbers deserve to be seen — reporting only the second quietly discards the runs closest to the null, the ones where stepwise selection picked nothing, and pushes the figure a few points further in the “stepwise regression manufactures findings” direction. Incidentally, even on pure noise, 92% of runs left stepAIC holding at least one variable.

Why this happens, and why it matters so much:

The selection is itself a multiple comparison, and the output does not account for it. The p-values in the final model are computed as though that set of variables had been specified in advance. They were not — they are the winners of a contest among many candidates. Every p-value, confidence interval and R² is therefore over-optimistic.
The coefficients are inflated by selection bias. The variables that survive are the ones whose effect happened to come out large in this sample, so the surviving coefficients sit systematically away from 0. The literature calls this Freedman’s paradox.
The result is extremely unstable. Refit on a different random subsample and the selected set is often very different. When two papers using the same method report different “independent risk factors”, the difference is frequently nothing but sampling.
It answers a question that is usually not yours. Stepwise selection is prediction-oriented variable screening, whereas most clinical papers want to estimate the effect of some exposure — a causal question, where the variables have to be chosen from a causal diagram.

Worth noting from the first simulation: smoking, a variable that genuinely exists, was retained every time (100%), with a mean coefficient of -355.5 against -360.7 in the full model. Stepwise regression is not incapable of finding real things; the problem is that it picks up a pile of false ones at the same time, and the output cannot tell them apart.

So how should the variables be chosen

If your question is	Do this
Estimating the effect of an exposure (most clinical research)	Let a causal diagram (DAG) decide which variables to adjust for, and write it into the protocol in advance. Do not let p-values decide what stays.
Building a prediction model	Work out the required sample size first, then use penalised methods (LASSO, ridge, elastic net) with internal validation (bootstrap or cross-validation) and external validation. Report following TRIPOD.
Exploratory work with no prior hypothesis	Legitimate, but label the whole analysis exploratory, do not report p-values as conclusions, and expect to need a second dataset to confirm anything.

All three share one feature: whether a variable stays is not decided by a p-value from the same dataset.

Run it yourself

library(MASS)
data(birthwt, package = "MASS")
bw <- birthwt
bw$race_f  <- factor(bw$race,  levels = 1:3, labels = c("White", "Black", "Other"))
bw$smoke_f <- factor(bw$smoke, levels = 0:1, labels = c("No", "Yes"))

fit <- lm(bwt ~ age + lwt + smoke_f + race_f + ht + ui, data = bw)

par(mfrow = c(2, 2)); plot(fit)          # the four standard diagnostic plots

# ⚠️ Those four plot against the fitted value; curvature in one variable gets diluted
plot(bw$age, residuals(fit)); lines(lowess(bw$age, residuals(fit)))

# VIF, computed from the definition (no need for car)
X <- model.matrix(fit)[, -1]
sapply(seq_len(ncol(X)), function(j) 1 / (1 - summary(lm(X[, j] ~ X[, -j]))$r.squared))

# Collinearity demo: add the same weight again, converted to kilograms
# ⚠️ The conversion needs a little error, or the two columns are exactly proportional,
#    R calls them aliased, the coefficient is NA and no VIF can be computed
bw$lwt_kg <- bw$lwt * 0.4536 + rnorm(nrow(bw), 0, 0.5)
fit_col <- lm(bwt ~ age + lwt + lwt_kg + smoke_f + race_f + ht + ui, data = bw)
summary(fit_col)$coefficients[c("lwt", "lwt_kg"), ]

# Influential points: find them, then actually refit and see whether the conclusion moves
cd <- cooks.distance(fit)
worst <- which.max(cd)
coef(fit)["smoke_fYes"]
coef(lm(formula(fit), data = bw[-worst, ]))["smoke_fYes"]

# What stepwise selection picks up: throw in 10 random variables and try it once
d <- bw[, c("bwt", "age", "lwt", "smoke_f", "race_f", "ht", "ui")]
for (j in 1:10) d[[paste0("z", j)]] <- rnorm(nrow(d))
summary(stepAIC(lm(bwt ~ ., data = d), direction = "both", trace = 0))

Verified with R 4.6.0 and MASS 7.3.65. The car package is not installed, so the VIFs on this page are computed from the definition; car::vif() returns GVIF for factors, which is a different number.

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.stats.outliers_influence import variance_inflation_factor

bw = sm.datasets.get_rdataset("birthwt", "MASS").data
bw["race_f"]  = bw["race"].map({1: "White", 2: "Black", 3: "Other"})
bw["smoke_f"] = bw["smoke"].map({0: "No", 1: "Yes"})

fit = smf.ols("bwt ~ age + lwt + C(smoke_f) + C(race_f) + ht + ui", data=bw).fit()

infl = fit.get_influence()
cooks = infl.cooks_distance[0]
lev   = infl.hat_matrix_diag

X = np.asarray(fit.model.exog)
vif = [variance_inflation_factor(X, i) for i in range(1, X.shape[1])]

# statsmodels has no built-in stepwise — which is a good thing

statsmodels' OLSInfluence provides cooks_distance / hat_matrix_diag / resid_studentized; VIF lives in statsmodels.stats.outliers_influence.

Common misuses

Misuse	Why it is wrong
Fitting a regression and looking at no diagnostics at all	The p-values and confidence intervals are only as good as the assumptions
Declaring the specification correct on the strength of the four standard plots	Curvature in one variable is diluted inside the fitted value; plot residuals variable by variable
Using a normality test to decide whether linear regression is allowed	The assumption is about the residuals, not the outcome, and it matters least in large samples
Dropping variables because the VIF is high	Collinearity only harms the variables caught up in it; if yours is not among them, leave it alone
Treating VIF 5 or 10 as a hard threshold	Both are conventions; the question is whether the standard error still answers your question
Comparing a car::vif() GVIF with a column-by-column VIF	They are defined differently and the values are not comparable
Deleting a high-leverage point as an outlier	High leverage is not influence; look at Cook’s distance
Cleaning data automatically with a Cook’s D threshold	Check whether it is a data error instead, and run a sensitivity analysis
Deleting influential points without declaring it	Data manipulation
Carrying the linear model’s residual plots over to logistic regression	Raw residuals from a 0/1 outcome only ever line up in two bands
Claiming good calibration because Hosmer-Lemeshow was not significant	That means “no evidence it is broken”; both the grouping and the sample size drive the result
Selecting variables by stepwise regression and reporting the p-values as usual	The selection is a multiple comparison; the p-values, CIs and R² are all over-optimistic
Calling the output of stepwise regression “independent risk factors”	That set of variables changes with the subsample
Screening by univariable p-value before building the multivariable model	The same problem in another form, and it invalidates the inference just as thoroughly

Reproducing every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B2-04-model-diagnostics.R

Read the figure

The answer comes from the same statistical output that produced this page's figures, not from a number typed in beside them.

Maternal weight is entered twice in the same model, once in pounds and once in kilograms. What happens to the lwt term?

Show the answer and why

Correct answer: Its VIF becomes 771.1 - two near-identical variables explain each other

lwt's VIF goes from 1.2 to 771.1. "Same information, different units" is true of the data and false of the estimates: two near-identical variables explain each other, the model cannot decide which owns the effect, and the individual standard errors blow up. 771.7 is the model's largest VIF, sitting on the kilogram version rather than on this row. Note that R-squared and the residual standard deviation barely move - collinearity does not make the model predict worse, it makes individual coefficients unreadable.

The duplicated weight variable wrecks lwt's VIF. What happens to the smoking coefficient?

Show the answer and why

Correct answer: It becomes -358.5 - a shift of well under a percent, because it is uncorrelated with the pair

Smoking moves from -360.7 to -358.5, a shift of well under one percent. Collinearity touches only the terms that are highly correlated with each other; coefficients unrelated to them are unaffected - so a huge VIF is not a reason to rebuild the whole model, the question is whether the term you care about got caught up in it. 647.4 is the residual standard deviation rather than any coefficient, and it barely moves either, which is the same point about prediction.

This model's overall residual plot, normality test and homoscedasticity test all look clean. So what is regressing the residuals on each covariate's square meant to reveal?

Show the answer and why

Correct answer: Whether age is linear: p = 0.017, evidence against linearity

Passing the overall diagnostics only means the residuals carry no structure visible at a glance; splitting by variable is what exposes the significant quadratic term on age (0.017), meaning age may not enter linearly. 0.749 belongs to maternal weight, which really is fine, and 0.632 is Tukey's non-additivity test, which asks about the model as a whole and misses age just as the other overall checks did. "All the diagnostics passed" often means only that nothing was split finely enough.

Chapters that use this method

Watch next

繁中臺大開放式課程 NTU OCW· 22 minIn Mandarin; lists the assumptions of linear regression one by one. Those are exactly what this page checks, so start by knowing what there is to check.

繁中NYCU OCW· 143 minNo verified video on this site is devoted to regression diagnostics; the residual and model-testing segments of this long lecture are the closest supplement.

Sources and licences

This page is original writing