Regression Diagnostics and Assumptions

A residual is observed minus predicted. Its pattern can reveal representation failure, variance changes, dependence, outliers, and data defects.

None of those messages reaches you through the coefficient table. A fitted slope, an R² and a standard error are summaries. Summaries are precisely the quantities that datasets with nothing in common can share. So the evidence lives in the individual residuals, not in the numbers printed above them. Hence the loop below. Fit on the designated training data. Predict out of fold. Stratify the residuals by prediction, feature range, time, entity and operational slice. Connect the visible pattern to a plausible mechanism. Then change one thing and retest.

Curvature against fitted values suggests missing nonlinear structure. A fan shape suggests error variance changes with scale. Long runs over time suggest dependence or drift, while isolated extreme residuals may indicate rare cases, data errors, or heavy-tailed noise.

A plot is evidence of mismatch, not a unique diagnosis. The converse is the sharper danger. The summary statistics that would let you skip the plot can agree exactly across shapes that share nothing.

The proof of that is a dinosaur. In 2017 Justin Matejka and George Fitzmaurice, at Autodesk Research, used simulated annealing to build the “Datasaurus Dozen”: twelve datasets that agree to two decimal places on an x-mean of 54.26, a y-mean of 47.83, an sd_x of 16.76, an sd_y of 26.93, and a Pearson r of -0.06. One of them plots as a dinosaur. The others plot as stars, lines, circles and slanted bands. Their abstract opens by naming the genre they were joining: “Datasets which are identical over a number of statistical properties, yet produce dissimilar graphs, are frequently used to illustrate the importance of graphical representations when exploring data.” What they added was a method. The CRAN package datasauRus now ships the twelve as a data frame of 1,846 rows and 3 variables, and its reference manual states the result in a line: “These 12 datasets are equal in standard measures: mean, standard deviation, and Pearson's correlation”.

Two consequences follow for residual work. Five agreeing numbers do not pin down a shape, so a diagnostic that checks only summaries will pass anything. And one visible shape does not pin down a mechanism. The plot tells you that a test is needed, not which test will pass.

The phrase “linear-regression assumptions” is incomplete unless the intended claim is stated. Point prediction needs the fitted relationship to survive the deployment distribution. It tolerates a misspecified model that is nonetheless useful, and it requires nothing about normal residuals. Classical uncertainty needs a defensible sampling and error model before an interval or a test means anything. That is where dependence and heteroscedasticity bite, and where software defaults stop counting as evidence. Causal interpretation needs an identification strategy connecting a coefficient to an intervention: control of confounding or randomization, correct temporal ordering, no inappropriate conditioning. The third column is the one that fails most quietly, because the failure prints as an ordinary table.

That failure has a name. Daniel Westreich and Sander Greenland called it the Table 2 fallacy in 2013. Working with causal diagrams, they showed that a single adjusted regression does not produce one kind of quantity. The adjustment set chosen to identify the exposure's effect generally does not identify the effects of the confounders. So direct-effect estimates get confused with total-effect estimates, and a covariate coefficient can remain confounded even when the exposure coefficient is not. The layout is what does the damage. “By presenting adjusted effect estimates for secondary risk factors alongside the adjusted effect estimate for the primary exposure, Table 2 suggests implicitly that all of these estimates can be interpreted similarly, if not identically.”

It is not a rare slip. A 2021 scoping review screened 1,358 articles from four oral-health journals published 2013-2018. Its results section reports the count: “The remaining 421 articles were eligible for full text reviewed of which, 189 (45%) committed Table 2 fallacy.” Nearly half the eligible studies in one field read every row of an adjusted model as an effect.

Each finding should suggest a testable next step rather than a decorative chart. The test is what separates a diagnosis from a story. The same picture admits several mechanisms, and only a refit distinguishes them.

Leverage describes how unusual a row is in feature space relative to the fitted design. A residual describes outcome mismatch. Influence describes how much the fitted model changes when a row is removed or downweighted. Two papers from 1977 turned each of those sentences into arithmetic you can compute and check.

The first one fixed the yardstick. Leverage is the diagonal element h_ii of the hat matrix, and because trace(H) = rank(X) = p, those diagonals have no freedom about their average. David C. Hoaglin and Roy E. Welsch wrote it in one line in January 1977: “The average size of a diagonal element of the hat matrix, then, is p/n .” The rule they built on top of it is still the one in use. Treat a point as high-leverage when h_ii exceeds 2p/n, twice the mean. Statistics Canada's Survey Methodology applies it word for word, calling an individual leverage “usually considered extreme if it is more than twice the mean”.

The other half arrived that February. R. Dennis Cook asked how far the least-squares estimate moves when observation i is deleted, measured in units of the normal-theory confidence ellipsoid for beta. His abstract states the construction: “A new measure based on confidence ellipsoids is developed for judging the contribution of each data point to the determination of the least squares estimate of the parameter vector in full rank linear regression models. It is shown that the measure combines information from the studentized residuals and the variances of the residuals and predicted values.” The measure is now universally called Cook's distance. It factors into exactly three quantities: the number of parameters p, the ith squared studentized residual, and the ratio V(fitted_i)/V(residual_i). D_i = (t_i² / p) × V(fitted_i)/V(residual_i). Influence is residual size and leverage multiplied, never either one alone. That is why the two dissociate so often. Statistics Canada's Survey Methodology puts the dissociation plainly: “some observations have high leverages and small residuals, while others have low leverages and large residuals”.

Read the algebra back into practice. A row with h_ii above 2p/n and a small studentized residual is well predicted and still anchoring the slope. A large residual at an ordinary feature location can hurt accuracy while barely moving the coefficients. Neither column alone tells you which rows your conclusions are resting on. The product does.

Run diagnostics on held-out or out-of-fold predictions whenever possible. Verify the rows first — units, impossible values, duplicates, target construction. Plot residuals against predictions for curvature, scale changes and boundary effects. Slice by group, time, missingness and operational regime. Measure sensitivity by refitting without the influential cases. Test one intervention on untouched evidence. Document what remains unresolved, and where the model should abstain or be monitored.

Step one and step four look like clerical work until you see what they catch. In April 2013 three economists obtained the working spreadsheet behind a widely cited result on public debt and growth. An averaging formula's range silently stopped short of the first five countries alphabetically: Australia, Austria, Belgium, Canada and Denmark. That coding error alone removed 25 Belgian country-years. Separate exclusions of available years removed more — 5 for Australia, 5 for Canada, 4 for New Zealand. The above-90%-debt bucket held 71 country-years where the corrected data contain 110. The excluded New Zealand years 1946-49 had grown 7.7%, 11.9%, -9.9% and 10.8%; with those gone, New Zealand entered the bucket through a single year, 1951, at -7.6%. Corrected, average real GDP growth above 90% debt/GDP is 2.2%. The published figure was -0.1%.

The ending is the part worth copying. On 5 May 2013 Carmen M. Reinhart and Kenneth S. Rogoff issued an errata that conceded the coding error and confirmed the postwar mean for that bucket rises by 0.3 points. It opens: “The authors wish to thank Thomas Herndon, Michael Ash, and Robert Pollin (HAP) for the careful attention and for pointing out an important correction to our original Growth in a Time of Debt.” No residual plot would have found this. Reading the rows found it. A sensitivity refit would have exposed it earlier: a headline that travels from -0.1% to 2.2% when a few dozen country-years come back was never resting on the model.

Heteroscedasticity-consistent standard errors can improve some inferential calculations when variance changes. They do not remove leakage. Nor do they remove unmodeled nonlinearity, dependent sampling, or distribution shift.

The cleanest demonstration is one where the mean function was correct by construction and the inference broke anyway. In the Quarterly Journal of Economics in 2004, three economists took roughly 900,000 CPS observations on women aged 25-50 for 1979-1999 — 50 states x 21 years, so 1,050 state-year cells — and invented placebo laws that affected nothing. Every "effect" recovered was therefore false by definition. Conventional standard errors recovered them at an extraordinary rate: “These conventional DD standard errors severely understate the standard deviation of the estimators: we find an "effect" significant at the 5 percent level for up to 45 percent of the placebo interventions.” Clustering the errors at the state-year level still left 44%. Their Monte Carlo replication produced 49%. Cameron and Miller reproduced the same Table 8 rejection rates independently in 2015: 0.063, 0.058, 0.080 and 0.115 for G = 50, 20, 10 and 6. Dependence among the residuals, not the fitted mean, was the entire problem.

Likewise, a high R² does not prove residuals are operationally acceptable, so inspect errors in the units and regions that drive decisions.

The corrections have a lineage worth knowing, and Achim Zeileis set it out in 2004. HC0 entered econometrics through White (1980), building on Eicker (1963). HC1, HC2 and HC3 “were suggested by MacKinnon and White (1985) to improve the performance in small samples”. Afterwards, “a more extensive study of small sample behaviour was carried out by Long and Ervin (2000) which arrive at the conclusion that HC3 provides the best performance in small samples as it gives less weight to influential observations”. R's sandwich package therefore defaults vcovHC to type = "HC3". Python's statsmodels ships the same estimator, documented as “MacKinnon and White’s (1985) heteroskedasticity robust standard errors” and defined as diag(e_i²/(1−h_ii)²) — the same h_ii that carries the leverage rule. Each variant names a different assumption about the variance. None of them names an assumption about the mean function, and none of them names an assumption about dependence.