Heteroscedasticity: how to detect it and fix it

Heteroscedasticity means the residuals of a regression are more spread out in some places than in others. It rarely biases the fitted line, but it makes the usual standard errors, confidence intervals and p-values unreliable, so the useful questions are how to see it and which fix matches its cause.

This post shows the residual plots to look at, explains what a formal test can and cannot tell you, and works one example in full on seeded synthetic data where the spread is known to grow with the predictor. The same rows are fitted twice through a deterministic engine: once on the original scale, and once after a log transform. A third fix, robust standard errors, is shown as code. Unequal spread is one of several assumptions behind a linear model; the overview of parametric test assumptions puts it in context.

What heteroscedasticity is, and what it is not

A linear regression models the mean of the outcome as a straight-line function of the predictors and treats everything else as noise. Ordinary least squares, the default fitting method, assumes that noise has the same variance at every value of the predictors. That is homoscedasticity. When the variance changes, typically growing as the fitted value grows, the residuals are heteroscedastic.

The usual cause is the measurement itself. Concentrations, incomes, reaction times, counts and anything else bounded at zero tend to vary more when they are large: the same proportional wobble is small in absolute terms around a small value and large around a large one. Other causes are a missing predictor whose effect differs across the range, a mixture of subgroups with different variability, and an outcome that is a count or a proportion modelled as if it were a continuous measurement.

Heteroscedasticity is a statement about the residuals, not the raw outcome. A skewed outcome can give perfectly even residuals once the predictors are accounted for, and a symmetric outcome can give a fan. It is also not the same as non-normal residuals, although the two often travel together, as the example below shows.

The same idea appears in group comparisons under another name. Unequal variances between groups in a t-test or ANOVA are heteroscedasticity with a categorical predictor, and Welch’s versions of those tests are the standard remedy, discussed in the guides on t-test vs Mann-Whitney and one-way ANOVA, Welch and Kruskal-Wallis.

What it does and does not damage

The coefficients first. If the model for the mean is correct, ordinary least squares estimates remain unbiased when the spread is unequal. Long and Ervin put it plainly in 2000: the estimates are unbiased, but the usual tests of significance are generally inappropriate.

The damage is to everything built from the standard errors. The usual formula assumes one common variance, so when the spread differs, the standard errors can be too small or too large, depending on where the noisy observations sit relative to the predictor’s range. Confidence intervals then cover the true value more or less often than their stated 95%, and p-values are off by the same mechanism. Ordinary least squares also stops being the most efficient estimator, because it gives the noisy observations the same weight as the precise ones.

Two consequences are easy to miss. Prediction intervals built on one common residual spread are too wide where the data are tight and too narrow where they are noisy, which matters if the model is used to predict. And unequal spread is often a symptom of a wrong mean model: the same multiplicative process that fans the residuals out also bends the relationship into a curve. In that case the coefficients of a straight line are not the problem to fix; the model is.

How to detect it

The residual plot comes first. Plot the residuals against the fitted values. With constant spread you see a horizontal band of roughly even width. With heteroscedasticity you usually see a funnel: a narrow band at one end and a wide one at the other. A curved band at the same time suggests the mean model is wrong as well.

fit <- lm(concentration ~ exposure, data = data)
plot(fit, which = 1)  # residuals vs fitted: look for a funnel or a curve
plot(fit, which = 3)  # scale-location: a rising trend means spread grows with the fitted value

The scale-location plot puts the square root of the absolute standardised residuals on the vertical axis, so growing spread shows as an upward trend, which is easier to judge than the width of a band. With several predictors, plot the residuals against each predictor as well, because the spread can depend on one predictor that the fitted values average away.

A formal test can back up the plot. The Breusch-Pagan test (Breusch and Pagan 1979) regresses the squared residuals on the predictors and asks whether they explain any of the variation. White’s 1980 paper gives a related direct test, which compares the robust and the ordinary estimates of the coefficients’ covariance. In R:

lmtest::bptest(fit)  # studentized Breusch-Pagan test by default

The Goldfeld-Quandt test (Goldfeld and Quandt 1965) suits spread that rises with one predictor: it orders the rows by that predictor, fits the regression separately to the low and high ends, and compares their residual variances with an F test (lmtest::gqtest() in R). In Python, statsmodels has het_breuschpagan(), studentized by default like bptest(), and fit(cov_type="HC3") for the robust errors in Fix 2. At small samples the Breusch-Pagan test can miss spread differences large enough to distort the standard errors, and at large samples it flags differences too small to matter. A non-significant result means unequal spread was not detected, not that the spread is equal. Long and Ervin went further and recommended that the decision to use robust standard errors should not be made by a screening test at all: if there is reason to suspect unequal spread, correct for it.

Worked example: spread that grows with the predictor

Real data never tell you for certain whether the spread is constant, so this example uses seeded synthetic data where the answer is known. One hundred observations of a predictor called exposure were drawn evenly between 0 and 10. The outcome, concentration, was generated with noise that is a constant percentage of its mean (log-normal, multiplicative error). That is the textbook process that produces a funnel: the larger the concentration, the larger its absolute wobble. The generating code, with the seed the engine run used:

RNGkind("Mersenne-Twister", "Inversion", "Rejection")
set.seed(20261022)
n <- 100L
exposure <- round(stats::runif(n, min = 0, max = 10), 2)
log_concentration <- 2 + 0.1 * exposure + stats::rnorm(n, sd = 0.35)
concentration <- round(exp(log_concentration), 3)
data <- data.frame(exposure, concentration)  # the name the other code blocks use

Seeded and run in R; a linear regression of concentration on exposure was chosen before the run.

The result on the original scale

Quantity	Value
Observations	100
Slope for exposure	1.488 per unit
95% confidence interval for the slope	1.085 to 1.891
Standard error of the slope (ordinary)	0.2032
t statistic	7.322 on 98 degrees of freedom
p-value	< 0.001
R-squared	0.3536
Shapiro-Wilk on the residuals	W = 0.9199, p < 0.001

Each unit of exposure was associated with an average increase of 1.488 in concentration, with a 95% confidence interval from 1.085 to 1.891, and the straight line explained about 35% of the variation (R-squared 0.3536). Strictly, the line is itself an approximation: with multiplicative noise the true mean of concentration is an exponential curve in exposure, mild over this range, so 1.488 describes the best-fitting line, not a parameter of the process.

The catch is that the interval is computed from the ordinary standard error, 0.2032, which assumes one common residual variance. These data were built to break that assumption. Rerun the generating code and draw the residual plot, and the residuals fan out as the fitted values grow. So the interval above rests on an assumption we know is false, and its stated 95% coverage cannot be taken at face value. That does not mean it is too narrow; as Fix 2 explains, the error can go either way.

What the engine’s checks said, and what they did not

The engine’s linear regression runs a normality check on the residuals, not a spread check. Here it failed: Shapiro-Wilk gave W = 0.9199 with p < 0.001, and the engine raised a non-normality flag (its wording refers to the outcome, but the check ran on the residuals) that suggests considering a transformation of the outcome if the skew is strong. The engine also noted that 5 of 100 observations had a Cook’s distance above the conventional 4/n screening line, the largest 0.194.

Neither flag names heteroscedasticity, but both are its fingerprints here. Multiplicative noise produces residuals that are skewed and fan-shaped at once, and the observations at the noisy end have the largest residuals, so they are the ones that look influential. The lesson is that a failed normality check in a regression is worth reading as “look at the residual plot”, not as “switch to a different test”. The engine did not test the spread itself, and a normality check cannot stand in for that plot.

Fix 1: transform the outcome

When the noise is proportional to the level, the logarithm turns a multiplicative process into an additive one, and the spread becomes even. The same 100 rows, fitted with log(concentration) as the outcome:

Quantity	Original scale	Log scale
Slope for exposure	1.488	0.1182
95% confidence interval	1.085 to 1.891	0.09117 to 0.1453
Standard error of the slope	0.2032	0.01363
t statistic	7.322	8.672
R-squared	0.3536	0.4342
Shapiro-Wilk on the residuals	W = 0.9199, p < 0.001	W = 0.9925, p = 0.8587
Largest Cook’s distance	0.194	0.078

On the log scale the residual normality check did not detect a departure (W = 0.9925, p = 0.8587), the largest Cook’s distance fell from 0.194 to 0.078 (5 rows still crossed the 4/n line), and the line explained more of the variation (R-squared 0.4342 against 0.3536). A residual plot of the log fit shows a band of roughly even width. Because these data were generated on the log scale, we also know the log fit is the right model, and its interval, 0.09117 to 0.1453, contains the slope used to generate the data.

The coefficient changes meaning. On the original scale, 1.488 is an additive change in concentration per unit of exposure. On the natural-log scale, 0.1182 is a change in log concentration, which means a multiplicative change in concentration itself: compute exp(0.1182) in R for the factor per unit, and exponentiate both interval limits the same way.

Three cautions. A transform changes the question as well as the spread: the log model describes the median or geometric mean, so back-transformed predictions are not arithmetic means. The R-squared comparison across scales is only a rough guide, because the two models explain variation in different outcomes. And a log transform suits noise that grows with the mean; it is the wrong fix when the mean is genuinely linear and only the spread varies, because it bends a straight relationship into a curve.

Fix 2: heteroscedasticity-consistent standard errors

When the straight-line mean is right and only the spread is unequal, keep the coefficients and replace the standard errors. White (1980) showed how to estimate the covariance of the coefficients consistently without modelling how the variance changes. Long and Ervin’s simulations found that the asymptotic version most software used (HC0) often gives incorrect inferences when there are 250 observations or fewer, while small-sample versions, especially one called HC3, work well even at 25; they recommended HC3 at 250 or fewer. The R package sandwich implements it (Zeileis 2004); the code below needs the sandwich and lmtest packages installed:

fit <- lm(concentration ~ exposure, data = data)
vc  <- sandwich::vcovHC(fit, type = "HC3")
lmtest::coeftest(fit, vcov. = vc)  # same estimates, robust standard errors, t and p
lmtest::coefci(fit, vcov. = vc)    # robust 95% confidence intervals

Two facts about this fix. First, it never changes a coefficient: the slope stays at exactly the ordinary least squares value, and only its standard error, interval and p-value change. That makes it the right tool when the estimate is meaningful on the original scale. Second, robust standard errors are not always larger than the ordinary ones. Whether they come out larger or smaller depends on how the unequal spread lines up with the predictor (in a simple regression, roughly whether the noisiest rows sit at the extremes of the predictor or near its middle), so run the code on your data rather than assume the ordinary interval is too narrow.

Fix 3: model the variance

When the cause of the unequal spread is known, model it rather than work around it. Weighted least squares gives each observation a weight inversely proportional to its variance, which restores efficiency, but only if the weights are right; guessed weights trade one assumption for another. A model built for the outcome type is often the cleaner answer: counts whose variance grows with the mean belong in a Poisson or negative binomial model (see Poisson vs negative binomial regression), and positive continuous outcomes with a constant coefficient of variation suit a gamma model with a log link.

Which fix to use

• Spread grows with the level, and the outcome is positive: try a log transform first and check the residual plot again. It often fixes the curve and the spread together; in the example it evened out the spread.
• The mean is linear and meaningful on the original scale, and only the spread varies: keep the model and report HC3 standard errors, especially at 250 observations or fewer.
• The outcome is a count, a proportion or a strictly positive skewed measurement: use a model whose variance function fits the outcome.
• In every case: decide from the design and the residual plot, not from whether a test crossed 0.05.

How to report it

Say what you checked, what you saw, and what you did. For the log-scale example:

A linear regression of log concentration on exposure was fitted to 100 observations. The residual-versus-fitted plot on the original scale showed spread increasing with the fitted value, so the outcome was log-transformed; on the log scale the spread was approximately constant. The slope was 0.1182 per unit of exposure on the natural-log scale (95% CI 0.09117 to 0.1453).

If you use robust standard errors instead, name the estimator (HC3) and the software, and report the robust interval rather than the ordinary one.

Rerunning the example

The engine output and the R script below the article come from the original-scale fit. The script reads the bundled data, fits the regression with lm() and prints its summary, which reproduces the slope, its standard error, the t statistic, the p-value and R-squared quoted above, and it recomputes the sensitivity calculation the engine reports. It does not print the confidence interval (call confint() on the lm() fit for that), rerun the Shapiro-Wilk or Cook’s distance checks, draw the residual plots, fit the log-scale model or compute robust standard errors: for those, run the code blocks in this post on the same data. The log-scale numbers come from a second engine run on the same seeded rows. The guide to reproducible statistical analysis explains what a reproducible analysis package should contain, and the linear regression method reference describes the method itself.

# Results

- Test: Linear regression
- n: 100
- statistic: 7.322
- p-value: < 0.001
- effect size: R-squared = 0.3536
- 95% CI: [1.085, 1.891]
- seed: 20261022

## Assumption checks
- normality (Shapiro-Wilk): passed=FALSE (evidence of non-normality)

## Power (sensitivity, not observed power)
- minimum detectable effect at n=100: Cohen's f2 = 0.08006
- to detect a conventional medium effect: 55 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/heteroscedasticity.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("concentration", "exposure")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- summary(lm(concentration ~ exposure, data = data))
print(result)

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