Linear regression assumptions and how to check them

Linear regression assumes that the mean of the outcome is a straight-line function of the predictors, that the errors are independent, that they have the same spread everywhere, and that they are roughly normal. Those four do not matter equally: a wrong mean structure biases the coefficients themselves, while unequal spread, dependence and non-normality mostly damage the standard errors, p-values and intervals.

This post takes each assumption in that order: what breaks when it fails, the R check, and the fix. It then works one example in full: a four-predictor regression on R’s LifeCycleSavings data, run through a deterministic engine, with its assumption output read in context.

The overview of parametric test assumptions covers the ideas shared with t-tests and ANOVA. Regression adds the shape of each relationship and predictors that overlap.

The model in one line

A linear regression says that each observed outcome equals a mean plus an error:

outcome = b0 + b1 x1 + b2 x2 + … + error

The part before the error is the mean structure. One family of assumptions is about it: is the right formula on the right-hand side? The other is about the errors: are they independent, equally spread and roughly normal? This split is what decides whether a violation biases your estimates or only your uncertainty about them.

Fitting the line by least squares needs no assumption. The assumptions come in when you want the coefficients to mean something and the standard errors, p-values and intervals to be right.

The assumptions at a glance

Assumption	If it fails	Main check	Usual remedy
Linearity (correct mean structure)	Coefficients are biased: they estimate the wrong quantity	Residuals vs fitted; component plus residual plots	Transform, add a curve (polynomial or spline), add the missing interaction
Independence of errors	Estimates usually unbiased; standard errors often too small	The study design; residuals in time or cluster order	Model the clustering (mixed model) or use cluster-robust errors
Constant error variance	Estimates unbiased; standard errors wrong in either direction	Residuals vs fitted; scale-location plot	Robust (HC3) standard errors; a transform; weighted least squares
Normal errors	Estimates unbiased; small-sample tests and prediction intervals affected	Normal Q-Q plot of residuals	Usually nothing in larger samples; bootstrap or a transform in small ones

Two things often listed as assumptions are not assumptions of the model, but you still check them:

• Multicollinearity. Overlap between predictors does not bias anything. It makes the affected coefficients imprecise. Only perfect collinearity breaks the model, because then a coefficient cannot be estimated at all.
• Influential observations. A few points can move the fitted line a long way. That is a property of your data, not a violated assumption, but it affects how far to trust a coefficient.

Under everything sits a design condition no residual plot can detect: the predictors should be measured without much error, and the model should contain the variables the question needs. Leave out an important confounder and the model estimates a real association that is not the effect you wanted. Econometrics calls the general version endogeneity: a predictor that is correlated with the error term, through an omitted variable, measurement error or the outcome feeding back on the predictor. It biases the coefficient, and the residuals cannot show it.

The outcome should be continuous (for counts, proportions or binary outcomes use the families linked below). Predictors can be continuous or categorical, with categories entered as indicator variables. No single minimum sample size makes the assumptions hold; what matters is enough observations per estimated coefficient for usefully narrow intervals, and checking the plots in the sample you have.

Assumption 1: linearity, the one that biases the estimates

Linearity means that the average outcome changes by the same amount for each one-unit step in a predictor, holding the others fixed. If the true relationship curves, levels off or depends on another predictor, a straight line fitted through it describes no real feature of the data. The slope you report is then an average over a shape the model does not have.

That is why linearity comes first: the other failures leave the coefficients right and the uncertainty wrong, while this one makes the coefficients wrong, and a larger sample only makes you more confident about the wrong number.

How to check it. Plot the residuals against the fitted values. With a correct mean structure the points scatter evenly around zero with no pattern; a U shape or an arch says the model is missing a curve. With several predictors, a component plus residual plot for each one shows whether that predictor’s relationship with the outcome is straight once the others are accounted for.

fit <- lm(sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings)
plot(fit, which = 1)    # residuals vs fitted
car::crPlots(fit)       # one component plus residual plot per predictor

What to do if it fails. Change the mean structure, not the test. A log transform of a predictor whose effect levels off, a quadratic term, a natural spline (splines::ns(x, df = 3)), or an interaction between two predictors whose effects depend on each other. If the outcome is a count or a proportion, the curvature may be telling you that a linear model is the wrong family; a Poisson or negative binomial model or a logistic regression fits those outcomes directly.

Assumption 2: independence of the errors

Independence means that knowing one observation’s error tells you nothing about another’s. It fails when observations come in clusters (patients within clinics, pupils within schools, repeated measurements on the same person) or in a sequence where neighbours resemble each other, such as daily or yearly series.

When errors are positively correlated within clusters, and a predictor varies mostly between clusters, the data hold less information than the row count suggests. The coefficients are usually still unbiased if the mean structure is right, but standard errors computed as if every row were independent can be seriously too small (Moulton 1990), so intervals come out too narrow and p-values too small.

How to check it. Mostly from the design: is any unit measured more than once, or grouped with units that would make it alike? For data in time order, plot the residuals against time and look for runs; car::durbinWatsonTest(fit) tests for correlation between neighbouring residuals, but it only means something when the row order is meaningful.

What to do if it fails. Model the dependence. For repeated measurements or clusters, a mixed model with a random effect for the cluster is the usual choice; the guide to repeated measures ANOVA versus mixed models covers when each fits. Cluster-robust standard errors are an alternative when the clustering is a nuisance rather than part of the question.

Assumption 3: constant variance of the errors

Constant variance (homoscedasticity) means the residuals have the same spread at every fitted value. A common failure is a funnel, where the spread grows with the prediction, as it often does for quantities that cannot go below zero.

Unequal variance does not bias ordinary least squares estimates. It makes the usual tests of significance generally inappropriate (Long and Ervin 2000), in either direction: the usual standard error assumes one error variance for every row, so if the larger errors sit at extreme predictor values, where they move the slope most, it comes out too small; if they sit near the middle, too large.

How to check it. The residuals vs fitted plot again, and the scale-location plot, which shows the square root of the absolute standardised residuals against the fitted values. A flat band is what you want; a rising or falling trend is a sign of unequal spread. The Breusch-Pagan test (lmtest::bptest(fit)) puts a p-value on it, with the usual limits: in small samples it misses real problems, and in large ones it flags differences too small to matter.

plot(fit, which = 3)                              # scale-location
lmtest::bptest(fit)                               # Breusch-Pagan test
lmtest::coeftest(fit, vcov = car::hccm(fit, type = "hc3"))  # HC3 standard errors

What to do if it fails. The simplest fix is heteroscedasticity-consistent standard errors, which keep the same coefficients and recompute their standard errors without assuming equal variance (White 1980). Long and Ervin’s simulations found that the HC3 version worked well even in samples as small as 25, and they recommend not letting a screening test decide whether to use it. A log transform of the outcome often fixes a funnel too, but it changes what the coefficients mean. Weighted least squares fits when you know how the variance changes.

Assumption 4: normality of the errors

The normality assumption is about the errors, which you see through the residuals. It is not about the predictors, and it is not about the raw outcome. A skewed predictor is fine, and the model does not require multivariate normality of the variables. A skewed outcome can be fine too, as long as the residuals from the fitted model are not badly skewed. Testing every variable for normality before a regression checks the wrong thing; a review of twelve clinical psychology journals found papers that treated normality of the variables themselves as a requirement (Ernst and Albers 2017).

Normality is also the least important of the four. The estimates do not depend on it, and in large samples the tests and intervals stay approximately valid for almost any error distribution (Lumley and colleagues 2002). It matters in small samples with strong skew or heavy tails, and for prediction intervals for single new observations at any sample size.

How to check it. A normal Q-Q plot of the residuals (plot(fit, which = 2)). Points close to the line, with some wobble at the ends, are what moderate samples look like. A Shapiro-Wilk test on the residuals gives a p-value, but read it with the same care as any assumption test: a non-significant result means no departure was detected, not that the errors are normal, and in large samples the test flags departures that are too small to affect anything.

What to do if it fails. In a large sample, usually nothing, once linearity and constant variance are in order: a curved mean structure or unequal spread often shows up as apparent non-normality, and fixing those fixes it. In a small sample, a bootstrap of the coefficients, or a transform chosen for the reasons above, are the usual routes.

Collinearity and influence: not assumptions, still checked

Collinearity. When two predictors carry much of the same information, the model struggles to split the credit between them. The variance inflation factor (VIF) is how many times larger a coefficient’s variance is than it would be if that predictor were unrelated to the others.

Rules of thumb (VIF above 5, or above 10) are conventions. O’Brien (2007) found that values far above them do not by themselves discount a regression’s results or call for dropping a predictor, and that the usual cures can create worse problems than the one they treat. Dropping a variable to lower a VIF changes the question the model answers and can introduce exactly the omitted-variable bias described above. The practical test: are the intervals for the coefficients you care about narrow enough for your purpose?

Influence. Cook’s distance combines how far a point sits from the others in the predictors (leverage) with how badly the model fits it (the residual) into one measure of how much the coefficients would change without it (Cook 1977). A common screening value is 4 divided by the sample size; like the VIF cut-offs it is a convention, not a law. A point above it is a point to look at, never one to delete automatically.

car::vif(fit)                               # variance inflation factors
cd <- cooks.distance(fit)
which(cd > 4 / nobs(fit))                   # rows above the 4/n screening value
plot(fit, which = 5)                        # residuals vs leverage

Correct influential points that are data-entry errors. For genuine ones, fit the model with and without them and report both; if the conclusion changes, that dependence is itself a finding.

Worked example: savings rates across 50 countries

R’s LifeCycleSavings data hold one row per country for 50 countries, with each variable averaged over 1960 to 1970: the aggregate personal savings ratio (sr), the percentage of the population under 15 (pop15) and over 75 (pop75), real per-capita disposable income (dpi) and its growth rate (ddpi). Belsley, Kuh and Welsch (1980) used this model as a worked example in their book on regression diagnostics. The research question for the run: how is the savings ratio related to age structure, income and income growth across countries?

The design settles two questions first. The outcome is continuous and each country contributes one row, so there are no repeated measurements; whether trading partners behave as independent units is a judgement no residual check can settle, and the analysis treats them as independent. Because the data are country averages, every coefficient describes countries, not individuals.

The engine fitted the model on all 50 rows (none were missing) with 45 residual degrees of freedom. The coefficient table:

Predictor	Estimate	95% CI	p
pop15	-0.4612	-0.7525 to -0.1699	0.0026
pop75	-1.691	-3.874 to 0.4910	0.126
dpi	-0.0003369	-0.002212 to 0.001538	0.719
ddpi	0.4097	0.0145 to 0.8049	0.0425

The headline coefficient is for pop15: t(45) = -3.189, p = 0.0026. Holding the other three predictors fixed, countries with a one-point larger share of under-15s had a savings ratio lower by 0.4612 points on average, with an interval from about 0.17 to 0.75 points lower. Together the four predictors give R-squared = 0.3385: they account for about a third of the variation in savings ratios across these countries. These are associations in observational data: they describe how the variables moved together across these countries, not what would happen to a country’s savings if its age structure changed.

The table holds four coefficient tests with no correction for multiplicity. pop15 is the headline because it was listed first, not because it won; ddpi, at p = 0.0425, sits near the 5% line and is weak evidence on its own, and comparing the two p-values does not rank the predictors.

Reading the assumption output

The engine runs its checks on the fitted model and reports them beside the result.

Normality of the residuals. Shapiro-Wilk on the residuals gave W = 0.987, p = 0.852. No departure from normality was detected. That is not evidence that the errors are normal, and the Q-Q plot from plot(fit, which = 2) is the better look at the shape. With 45 residual degrees of freedom and no sign of heavy skew, normality is the least of this model’s concerns.

Collinearity. The engine flagged pop15 and pop75 with variance inflation factors above 5, its screening value. A country with many children tends to have few people over 75, and that explains the table. The pop75 estimate is large and negative, but its interval runs from -3.874 to 0.4910, wide enough to include zero and also effects more than twice the estimate. The data cannot tell the two age shares apart well, so they cannot say much about pop75 once pop15 is in the model. Dropping one to lower the VIF would change the question and load its association onto the other. If the question is about age structure as a whole, test the two shares together with an F test comparing the model with and without both:

anova(update(fit, . ~ . - pop15 - pop75), fit)   # joint F test for the two age shares

Influence. Three of the 50 observations had a Cook’s distance above the 4/n screening value, the largest being 0.268. The engine reports the count, not which rows, so the next step is to list them with the code above, check them against their source, and refit without them. It also noted values beyond Tukey’s interquartile-range fence in sr (one) and ddpi (three); those flag single variables, and the influence check is what says whether they matter for the model.

What the engine does not check. It reports no test of linearity, constant variance or independence for this model. Those need the plots and design reasoning above. If the HC3 intervals disagree materially with the ordinary ones for a coefficient you care about, distrust the ordinary ones.

Sensitivity, not observed power

With 50 countries, four predictors and the conventional 5% significance level, the sensitivity calculation shows the smallest effect the model had 80% power to detect: Cohen’s f-squared = 0.2651 for the four predictors jointly. Detecting a conventional medium-sized effect with the same power would take 85 countries. This describes what the design could detect, not the power for the observed effects, and is no basis for sizing a future study.

Writing it up

A results paragraph could read (the bracket marks the refit this example did not run):

A linear regression of the aggregate savings ratio on the percentage of the population under 15 and over 75, per-capita disposable income and its growth rate (n = 50 countries) found that a larger share of under-15s was associated with a lower savings ratio (b = -0.4612, 95% CI -0.7525 to -0.1699, t(45) = -3.189, p = 0.0026), adjusting for the other predictors. The four predictors together gave R-squared = 0.3385. The two age-structure variables were collinear (variance inflation factors above 5), so the estimate for the share over 75 was imprecise (95% CI -3.874 to 0.4910). Residual normality was not rejected (Shapiro-Wilk p = 0.852). Three observations exceeded the 4/n Cook’s distance screening value; [state here whether the conclusions held when the model was refitted without them].

Which violations bias the estimates

A wrong mean structure, an omitted confounder and substantial measurement error in a predictor bias the coefficients, and more data does not fix them. The other failures mainly damage standard errors, small-sample inference and prediction intervals.

Choosing whether linear regression is the right model at all comes before any of this: the guide on how to choose a statistical test starts from the outcome type, and the linear regression method reference lists what the engine computes for this model.

Rerunning the example

The example’s reproduce.R reads the bundled data file, refits the same model with summary(lm(sr ~ pop15 + pop75 + dpi + ddpi, data = data)) and runs the same sensitivity calculation with pwr::pwr.f2.test(). Running it with the pinned package versions prints the estimates, t statistics, p-values, R-squared and minimum detectable effect quoted above; confint() on the same model gives the intervals. It does not compute the 85-country figure or rerun the engine’s own checks (Shapiro-Wilk, variance inflation factors, Cook’s distances, Tukey fences); the diagnostic code above computes those from the same model. The reproducibility guide explains what a seeded, version-pinned script can and cannot guarantee.

# Results

- Test: Linear regression
- n: 50
- statistic: -3.189
- p-value: 0.0026
- effect size: R-squared = 0.3385
- 95% CI: [-0.7525, -0.1699]
- seed: 20261022

## Assumption checks
- normality (Shapiro-Wilk): passed=TRUE (no evidence against normality)

## Power (sensitivity, not observed power)
- minimum detectable effect at n=50: Cohen's f2 = 0.2651
- to detect a conventional medium effect: 85 total to detect f2=0.15
- alpha = 0.05, target power = 0.8

# Celeus - reproduce.R  (run from this package folder; bundled data in data/)
# Test: Linear regression  |  engine 0.30.0  |  catalog c39-d3cd
# Research context: not provided for this run.
# Reproduces the recorded result deterministically.

set.seed(20261022)
.data_path <- "data/linear-regression-assumptions.csv"
if (!file.exists(.data_path)) stop("Run reproduce.R from the package folder (the bundled data is in data/).")
data <- readr::read_csv(.data_path, show_col_types = FALSE)  # needs the readr package
clean_cells <- function (df) 
{
    ctrl_collapse <- "[\r\n\t\v\f]+"
    zero_width_chars <- intToUtf8(c(8203L, 65279L), multiple = TRUE)
    uni_space_chars <- intToUtf8(c(160L, 8239L, 8199L), multiple = TRUE)
    clean_text <- function(x) {
        out <- x
        for (ch in zero_width_chars) out <- gsub(ch, "", out, 
            fixed = TRUE)
        out <- gsub(ctrl_collapse, " ", out)
        for (ch in uni_space_chars) out <- gsub(ch, " ", out, 
            fixed = TRUE)
        trimws(gsub(" +", " ", out))
    }
    changed_cols <- character(0)
    cells_cleaned <- 0L
    retyped <- character(0)
    for (j in seq_along(df)) {
        col <- df[[j]]
        if (!is.character(col)) 
            next
        cleaned <- clean_text(col)
        n_changed <- sum(!is.na(col) & !is.na(cleaned) & col != 
            cleaned, na.rm = TRUE)
        if (n_changed > 0L) {
            changed_cols <- c(changed_cols, names(df)[j])
            cells_cleaned <- cells_cleaned + n_changed
        }
        df[[j]] <- cleaned
        non_blank <- cleaned[!is.na(cleaned) & nzchar(cleaned)]
        if (length(non_blank) && all(grepl("^[+-]?(0|[1-9][0-9]*)([.][0-9]+)?$", 
            non_blank))) {
            parsed <- suppressWarnings(as.numeric(cleaned))
            was_blank <- is.na(cleaned) | !nzchar(cleaned)
            lost <- (is.na(parsed) | is.infinite(parsed)) & !was_blank
            if (!any(lost)) {
                finite <- parsed[!is.na(parsed)]
                all_whole <- all(finite == floor(finite))
                in_range <- all(abs(finite) <= .Machine$integer.max)
                df[[j]] <- if (all_whole && in_range) 
                  as.integer(parsed)
                else parsed
                retyped <- c(retyped, names(df)[j])
            }
        }
    }
    names(df) <- clean_text(names(df))
    if (length(changed_cols)) {
        attr(df, "ingest_clean_columns") <- changed_cols
        attr(df, "ingest_clean_cells_cleaned") <- cells_cleaned
    }
    if (length(retyped)) 
        attr(df, "ingest_clean_retyped") <- retyped
    df
}
data <- clean_cells(data)
# Complete-case analysis: drop rows missing any analysis value.
.acols <- c("sr", "pop15", "pop75", "dpi", "ddpi")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- summary(lm(sr ~ pop15 + pop75 + dpi + ddpi, data = data))
print(result)

# Sensitivity: minimum detectable effect at n=50 (power=0.8, alpha=0.05). NOT observed power.
power <- pwr::pwr.f2.test(u = 4, v = 45, sig.level = 0.05, power = 0.8)
print(power)