Breusch-Pagan test in R: what it tells you

The Breusch-Pagan test asks one question of a fitted regression: does the spread of the residuals change in a straight-line way with the predictors? It answers with a p-value, and that p-value tells you whether a change in spread was detected, not whether it is big enough to matter.

This post shows what the test computes, how to run and read it in R on the built-in cars data, which of its two versions to use, and a seeded demonstration of its main weakness: with enough rows it flags a change in spread that, we argue below, is too small to matter for most purposes. The constant-spread assumption itself, and what to do about each assumption that fails, is covered in the parametric assumptions guide. For the wider problem of spotting non-constant spread and handling it, see heteroscedasticity in regression.

What the test computes

Ordinary least squares regression assumes the errors have the same variance at every value of the predictors. When that holds, the usual standard errors, confidence intervals and p-values for the coefficients are correct. When it fails (heteroscedasticity), the coefficients are still unbiased, but their standard errors can be too small or too large.

Breusch and Pagan (1979) built a test for this from the residuals alone. The idea fits in three steps:

1. Fit the model and keep the residuals.
2. Square them. A squared residual is a rough, one-row estimate of the error variance at that point.
3. Regress the squared residuals on the predictors (the auxiliary regression). If the spread is constant, the predictors should explain almost none of the variation in the squared residuals.

The test statistic measures how much the auxiliary regression explains. Under constant spread it follows, approximately, a chi-square distribution with as many degrees of freedom as there are predictors in the auxiliary regression, not counting the intercept.

Two consequences follow from that construction. The test looks only for spread that changes with the variables you give it, so a pattern tied to something else goes unseen. And it looks mainly for a linear trend in the variance; a funnel that opens and then closes can slip past it.

White’s (1980) test is the usual alternative. It regresses the squared residuals on the predictors, their squares and their cross-products, so it can catch curved patterns in the spread. The price is more terms, which dilutes the test when you only worry about a simple trend. Choose between them by the pattern you suspect; one is not a stricter version of the other.

Two versions, and why the default matters

There are two versions in common use, and R’s two best-known functions run different ones.

• The original (score) version. Breusch and Pagan’s statistic, and the closely related score test of Cook and Weisberg (1983), scale the auxiliary regression by the variance a squared residual would have if the errors were normal. car::ncvTest() computes this version; by default it uses the fitted values as the variance predictor, which is the same thing when the model has one predictor.
• The studentized version. Koenker (1981) proposed studentizing the test: the normal-theory scaling is replaced by one estimated from the squared residuals themselves, and the statistic becomes n times the R-squared of the auxiliary regression. lmtest::bptest() computes this version unless you set studentize = FALSE, and by default it uses the model’s own predictors as the variance predictors.

The difference matters when the errors are not normal. The reasoning is short: the original statistic divides by the variance a squared residual would have under normal errors, and heavy-tailed errors make squared residuals more variable than that, so the statistic is inflated and the test can reject for reasons that have nothing to do with spread. The studentized version estimates that variance from the data instead, which is why it is the safer default. If you report a Breusch-Pagan test, say which version you ran.

Running it in R on the cars data

The cars data hold 50 measurements of stopping distance (feet) and speed (miles per hour), recorded in the 1920s. Regressing distance on speed is a common teaching example for a spread check.

The model first. Run through the deterministic engine, the regression gives:

• slope for speed: 3.932 feet per mile per hour, 95% CI 3.097 to 4.768
• t = 9.464 on 48 degrees of freedom, p < 0.001
• R-squared = 0.6511, n = 50

The engine’s assumption check on this model is a Shapiro-Wilk test of the residuals, and it does detect non-normality: W = 0.945, p = 0.0215. The engine also notes 2 of the 50 rows with a Cook’s distance above the usual screening threshold. Neither is a spread test, but the first is a reason to prefer the studentized version here: the normal-error scaling of the original is in doubt.

This page does not show a Breusch-Pagan result for the cars model. The engine’s regression output does not include the test, and every number on this page comes from engine output, so the test is given here as code for you to run, not as a result:

fit <- lm(dist ~ speed, data = cars)

# Residuals against fitted values: look for a funnel
plot(fitted(fit), resid(fit), xlab = "Fitted distance", ylab = "Residual")
abline(h = 0, lty = 2)

lmtest::bptest(fit)                     # studentized (Koenker), the default
lmtest::bptest(fit, studentize = FALSE) # original Breusch-Pagan
car::ncvTest(fit)                       # score test, original scaling

# White-style: squares (and cross-products, with several predictors) as variance predictors
lmtest::bptest(fit, ~ speed + I(speed^2), data = cars)

How to read the output

bptest() prints three numbers: BP, the test statistic; df, the number of predictors in the auxiliary regression (1 here, for speed); and the p-value from the chi-square distribution on those degrees of freedom.

Read them in this order:

1. The p-value says whether a linear change in spread was detected at your chosen level. Below 0.05: the data are hard to square with constant spread. Above it: a change in spread was not detected. That is not evidence that the spread is constant, and with 50 rows the test can miss changes that would matter.
2. The version. The heading says studentized Breusch-Pagan test for Koenker’s version. If the two versions disagree on your model, and the residuals look non-normal, trust the studentized one.
3. The size. The statistic is n times the auxiliary R-squared, so for the studentized version you can recover that R-squared by dividing BP by n. It is the share of the variation in squared residuals that the predictors explain, which is a rough measure of how strongly the spread tracks them.

Then look at the residual plot: the test gives one number for one pattern, the plot shows which pattern is there.

A seeded demonstration: what the test sees, and what it misses

To show the test’s behaviour on data where the truth is known, the same auxiliary regression was run through the engine on three seeded synthetic data sets. In each, x is uniform between 0 and 10, the true line is y = 5 + 2x, and the errors are normal with a standard deviation that rises with x. The engine regressed the squared residuals on x. The t test of that slope is the F form of the studentized test (F is t squared); bptest() reports the chi-square form of the same regression, and at these sample sizes the two forms give nearly the same p-value.

Data	Residual SD across the range	n	Auxiliary R-squared	Slope t	p
Strong	1 rising to 11	50	0.1039	2.359	0.0225
Mild	10 rising to 11	200	0.001977	0.6262	0.532
Mild	10 rising to 11	20000	0.001291	5.084	< 0.001

Three lessons sit in that table.

A large change in spread can be only just detected in a small sample. In the strong data the variance grows more than a hundredfold across the range of x, and with 50 rows the test gives p = 0.0225. Squared residuals are noisy, so the auxiliary regression needs many rows to see a trend in them.

A mild change is not detected at 200 rows. With the standard deviation rising by a tenth, the test gives p = 0.532. The change is real, because it was built into the data. A non-significant result is exactly what “not detected” looks like.

The same mild change is flagged at 20000 rows. The auxiliary regression estimates the average squared residual at about 101.3 when x is 0, rising by 1.948 per unit of x (95% CI 1.197 to 2.699). Across the whole range that is an increase in the error variance of roughly a fifth (a standard deviation up by about a tenth, as the data were built), and x explains about 0.13 percent of the variation in the squared residuals. The p-value is below 0.001 all the same. Whether a change that size matters is a judgement this demonstration does not measure: we did not compare conventional and robust standard errors on these data. Our argument is from how they were built: with the error standard deviation differing by only a tenth between the two ends of x, every row carries nearly the same noise, so conventional standard errors should stay close to robust ones. The next section shows the check for your own model. The test is doing its job: it answers “is the spread exactly constant?”, and with 20000 rows even a small departure is likely to be detected.

What to do with a significant result

A significant Breusch-Pagan test is a prompt to look, not an instruction to change the model. Two checks tell you whether it matters for your conclusions:

# Conventional versus heteroscedasticity-consistent (HC3) standard errors
lmtest::coeftest(fit)
lmtest::coeftest(fit, vcov. = car::hccm(fit, type = "hc3"))

If the robust standard errors are close to the conventional ones, the change in spread is not affecting your inference. If they differ, report the robust ones. White (1980) introduced these heteroscedasticity-consistent standard errors, and Long and Ervin (2000) found in simulations that the HC3 version works well even in small samples, while the original (HC0) version often gives wrong inferences at 250 rows or fewer.

Long and Ervin also recommend that the decision to use robust standard errors should not be made by a screening test for heteroscedasticity. That is the practical place of the Breusch-Pagan test: it describes the residuals, and it can confirm what a plot suggests, but whether to use robust standard errors is better settled by the design and the question than by its p-value. The same logic is behind using Welch’s t-test by default rather than choosing it with a variance test, and behind Welch’s ANOVA for groups with unequal spread.

Other remedies change the model rather than its standard errors: a transformation of the outcome (a log when the spread grows with the mean), weighted least squares when you can model the variance, or a model whose outcome distribution expects the spread to change. Each changes the question being answered, so choose it for the data, not for the test. If you are still deciding which model fits your outcome, start with the test chooser.

Reporting it

A sentence that gives the reader what they need:

Residual spread was checked with a residual-versus-fitted plot and the studentized Breusch-Pagan test (lmtest::bptest, speed as the variance predictor). We report heteroscedasticity-consistent (HC3) standard errors alongside the conventional ones.

In a results line, give the statistic, its degrees of freedom and the exact p-value, in the form BP(df) = statistic, p = value, then say what you did next. Name the version, name the variance predictors, and report what you did about it rather than the p-value alone. Python’s statsmodels also defaults to Koenker’s version (het_breuschpagan() with robust=True); for other software, check which version the default computes before you compare outputs.

Rerunning the example

The script in the engine package reproduces the regression of distance on speed (coefficients, standard errors, p-values, R-squared) and the attached sensitivity power calculation. It does not print the confidence interval (run confint(fit) for that), recompute the Shapiro-Wilk check or the Cook’s distance note, or run the Breusch-Pagan test or the seeded demonstration. The reproducible analysis guide explains what such a package contains, and the linear regression reference lists the engine’s outputs for this method.

This example is about the method, not about cars or braking: the 1920s data are public and small.

# Results

- Test: Linear regression
- n: 50
- statistic: 9.464
- p-value: < 0.001
- effect size: R-squared = 0.6511
- 95% CI: [3.097, 4.768]
- seed: 20261022

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

## Power (sensitivity, not observed power)
- minimum detectable effect at n=50: Cohen's f2 = 0.1635
- 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/breusch-pagan-test.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("dist", "speed")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- summary(lm(dist ~ speed, 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 = 1, v = 48, sig.level = 0.05, power = 0.8)
print(power)