How to Perform the Breusch-Pagan Test in R

Key Insights

The Breusch-Pagan test detects heteroscedasticity by regressing squared residuals on your predictors—a significant result (p < 0.05) indicates non-constant variance that violates OLS assumptions.
Use bptest() from the lmtest package for the standard test, or ncvTest() from car for an alternative implementation with slightly different defaults.
When heteroscedasticity is detected, don’t abandon your model—use robust standard errors via the sandwich package to get valid inference without changing your coefficient estimates.

Introduction to Heteroscedasticity

Heteroscedasticity occurs when the variance of regression residuals changes across the range of predictor values. This violates a core assumption of ordinary least squares (OLS) regression: that errors have constant variance (homoscedasticity).

The practical consequences are serious. Your coefficient estimates remain unbiased, but your standard errors become unreliable. This means hypothesis tests and confidence intervals can’t be trusted. You might declare a variable statistically significant when it isn’t, or miss genuine effects because your standard errors are inflated.

Let’s visualize the difference between homoscedastic and heteroscedastic residuals:

set.seed(42)

# Generate homoscedastic data (constant variance)
x <- seq(1, 100, length.out = 200)
y_homo <- 2 + 0.5 * x + rnorm(200, mean = 0, sd = 5)

# Generate heteroscedastic data (variance increases with x)
y_hetero <- 2 + 0.5 * x + rnorm(200, mean = 0, sd = 0.5 * x)

# Create comparison plots
par(mfrow = c(1, 2))

# Homoscedastic
plot(x, y_homo, main = "Homoscedastic Residuals",
     xlab = "X", ylab = "Y", pch = 16, col = "steelblue")
abline(lm(y_homo ~ x), col = "red", lwd = 2)

# Heteroscedastic
plot(x, y_hetero, main = "Heteroscedastic Residuals",
     xlab = "X", ylab = "Y", pch = 16, col = "steelblue")
abline(lm(y_hetero ~ x), col = "red", lwd = 2)

par(mfrow = c(1, 1))

In the homoscedastic plot, points scatter evenly around the regression line. In the heteroscedastic plot, the spread fans out as X increases—a classic pattern in real-world data involving income, prices, or counts.

What is the Breusch-Pagan Test?

The Breusch-Pagan test, developed by Trevor Breusch and Adrian Pagan in 1979, provides a formal statistical test for heteroscedasticity. The logic is elegant: if variance is truly constant, then squared residuals shouldn’t be predictable from your independent variables.

The test works by:

Fitting your original regression and extracting residuals
Squaring those residuals
Regressing squared residuals on the original predictors
Testing whether this auxiliary regression explains significant variance

Null hypothesis (H₀): Homoscedasticity—residual variance is constant across all values of the predictors.

Alternative hypothesis (H₁): Heteroscedasticity—residual variance depends on one or more predictors.

The test statistic follows a chi-squared distribution with degrees of freedom equal to the number of predictors in the auxiliary regression. A large test statistic (small p-value) suggests the squared residuals are systematically related to your predictors, indicating heteroscedasticity.

Setting Up Your R Environment

You need two packages for heteroscedasticity testing and remediation:

# Install packages if not already installed
install.packages("lmtest")
install.packages("car")
install.packages("sandwich")

# Load the packages
library(lmtest)
library(car)
library(sandwich)

The lmtest package provides bptest() for the Breusch-Pagan test. The car package offers ncvTest() as an alternative. The sandwich package supplies robust standard error estimators for when you detect problems.

Let’s prepare two datasets—one we know has heteroscedasticity and one that doesn’t:

set.seed(123)

# Dataset with heteroscedasticity
n <- 150
x1 <- runif(n, 1, 50)
x2 <- runif(n, 10, 100)

# Error variance increases with x1
hetero_error <- rnorm(n, mean = 0, sd = 2 * sqrt(x1))
y_hetero <- 10 + 2 * x1 + 0.5 * x2 + hetero_error

df_hetero <- data.frame(y = y_hetero, x1 = x1, x2 = x2)

# Dataset with homoscedasticity
homo_error <- rnorm(n, mean = 0, sd = 10)
y_homo <- 10 + 2 * x1 + 0.5 * x2 + homo_error

df_homo <- data.frame(y = y_homo, x1 = x1, x2 = x2)

We’ve created two scenarios with identical underlying relationships but different error structures. This lets us see how the test behaves under both conditions.

Running the Breusch-Pagan Test

The workflow is straightforward: fit your linear model first, then pass it to bptest().

# Fit linear models
model_hetero <- lm(y ~ x1 + x2, data = df_hetero)
model_homo <- lm(y ~ x1 + x2, data = df_homo)

# Run Breusch-Pagan test on heteroscedastic model
bp_hetero <- bptest(model_hetero)
print(bp_hetero)

Output:

	studentized Breusch-Pagan test

data:  model_hetero
BP = 18.432, df = 2, p-value = 9.94e-05

# Run Breusch-Pagan test on homoscedastic model
bp_homo <- bptest(model_homo)
print(bp_homo)

Output:

	studentized Breusch-Pagan test

data:  model_homo
BP = 1.287, df = 2, p-value = 0.5254

The bptest() function accepts several useful arguments:

# Test against specific variables only
bptest(model_hetero, ~ x1, data = df_hetero)

# Use original (non-studentized) version
bptest(model_hetero, studentize = FALSE)

# Test against fitted values instead of predictors
bptest(model_hetero, ~ fitted(model_hetero))

By default, bptest() uses the studentized version (Koenker’s variant), which is more robust to non-normality. The studentize = FALSE option gives you the original Breusch-Pagan formulation, which assumes normally distributed errors.

Interpreting the Results

Let’s break down what the output tells us:

# Extract components from the test result
bp_result <- bptest(model_hetero)

cat("Test Statistic (BP):", bp_result$statistic, "\n")
cat("Degrees of Freedom:", bp_result$parameter, "\n")
cat("P-value:", bp_result$p.value, "\n")

Test Statistic (BP): This is the chi-squared statistic. Larger values indicate stronger evidence of heteroscedasticity.

Degrees of Freedom: Equal to the number of predictors used in the auxiliary regression (2 in our case: x1 and x2).

P-value: The probability of observing this test statistic if the null hypothesis (homoscedasticity) were true.

Here’s how to make decisions programmatically:

interpret_bp_test <- function(model, alpha = 0.05) {
  bp <- bptest(model)
  
  cat("Breusch-Pagan Test Results\n")
  cat("==========================\n")
  cat("Test statistic:", round(bp$statistic, 4), "\n")
  cat("Degrees of freedom:", bp$parameter, "\n")
  cat("P-value:", format.pval(bp$p.value, digits = 4), "\n\n")
  
  if (bp$p.value < alpha) {
    cat("Decision: REJECT null hypothesis at alpha =", alpha, "\n")
    cat("Conclusion: Evidence of heteroscedasticity detected.\n")
    cat("Action: Consider robust standard errors or WLS.\n")
  } else {
    cat("Decision: FAIL TO REJECT null hypothesis at alpha =", alpha, "\n")
    cat("Conclusion: No significant evidence of heteroscedasticity.\n")
    cat("Action: OLS standard errors are likely valid.\n")
  }
  
  invisible(bp)
}

# Apply to both models
interpret_bp_test(model_hetero)
interpret_bp_test(model_homo)

For the heteroscedastic model, the p-value of 9.94e-05 is far below 0.05, so we reject the null hypothesis. For the homoscedastic model, the p-value of 0.525 provides no evidence against constant variance.

You can also use the ncvTest() function from the car package, which performs a similar test:

ncvTest(model_hetero)

The results will be comparable, though not identical due to slight implementation differences.

Addressing Heteroscedasticity

When you detect heteroscedasticity, you have several options. The most practical is using heteroscedasticity-consistent (HC) standard errors, often called “robust” standard errors.

# Original model summary with potentially invalid SEs
summary(model_hetero)

# Robust standard errors using sandwich estimator
robust_se <- coeftest(model_hetero, vcov = vcovHC(model_hetero, type = "HC3"))
print(robust_se)

The vcovHC() function computes heteroscedasticity-consistent covariance matrices. The “HC3” type is recommended for small to moderate samples—it’s more conservative than alternatives.

Compare the standard errors:

# Side-by-side comparison
original_se <- summary(model_hetero)$coefficients[, "Std. Error"]
robust_se_vals <- sqrt(diag(vcovHC(model_hetero, type = "HC3")))

comparison <- data.frame(
  Coefficient = names(coef(model_hetero)),
  Original_SE = round(original_se, 4),
  Robust_SE = round(robust_se_vals, 4),
  Ratio = round(robust_se_vals / original_se, 2)
)
print(comparison)

You’ll typically see robust standard errors differ from OLS standard errors when heteroscedasticity is present. The ratio tells you whether OLS was overestimating or underestimating uncertainty.

For weighted least squares, you need to model the variance structure:

# Estimate weights (inverse of variance)
# If variance proportional to x1:
weights <- 1 / fitted(lm(abs(residuals(model_hetero)) ~ x1, data = df_hetero))^2

# Fit weighted least squares
model_wls <- lm(y ~ x1 + x2, data = df_hetero, weights = weights)
summary(model_wls)

# Verify heteroscedasticity is reduced
bptest(model_wls)

WLS requires you to correctly specify how variance changes with predictors. If you get it wrong, you can make things worse. Robust standard errors are safer when you’re uncertain about the variance structure.

Conclusion

The Breusch-Pagan test is your first-line diagnostic for heteroscedasticity in linear regression. Run it routinely as part of model validation—it takes one line of code and can save you from drawing incorrect conclusions.

Key points to remember:

A significant result (p < 0.05) means heteroscedasticity is present. Your coefficients are still unbiased, but standard errors and tests are unreliable.
The studentized version (studentize = TRUE, the default) is more robust to non-normal errors.
When heteroscedasticity is detected, robust standard errors are usually the simplest fix.

The test has limitations. It can be sensitive to model misspecification and may have low power in small samples. Consider supplementing with visual diagnostics—plot residuals against fitted values and each predictor.

Alternative tests worth knowing: White’s test is more general and doesn’t assume a specific form of heteroscedasticity. The Goldfeld-Quandt test is useful when you suspect variance changes at a specific point in your data. Both are available in R through the lmtest package.