A Breusch-Pagan Test is used to determine if heteroscedasticity is present in a regression analysis.
This tutorial explains how to perform a Breusch-Pagan Test in R.
Example: Breusch-Pagan Test in R
In this example we will fit a regression model using the built-in R dataset mtcars and then perform a Breusch-Pagan Test using the bptest function from the lmtest library to determine if heteroscedasticity is present.
Step 1: Fit a regression model.
First, we will fit a regression model using mpg as the response variable and disp and hp as the two explanatory variables.
#load the dataset data(mtcars) #fit a regression model model <- lm(mpg~disp+hp, data=mtcars) #view model summary summary(model) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 30.735904 1.331566 23.083 < 2e-16 *** disp -0.030346 0.007405 -4.098 0.000306 *** hp -0.024840 0.013385 -1.856 0.073679 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 3.127 on 29 degrees of freedom Multiple R-squared: 0.7482, Adjusted R-squared: 0.7309 F-statistic: 43.09 on 2 and 29 DF, p-value: 2.062e-09
Step 2: Perform a Breusch-Pagan Test.
Next, we will perform a Breusch-Pagan Test to determine if heteroscedasticity is present.
#load lmtest library library(lmtest) #perform Breusch-Pagan Test bptest(model) studentized Breusch-Pagan test data: model BP = 4.0861, df = 2, p-value = 0.1296
The test statistic is 4.0861 and the corresponding p-value is 0.1296. Since the p-value is not less than 0.05, we fail to reject the null hypothesis. We do not have sufficient evidence to say that heteroscedasticity is present in the regression model.
What To Do Next
If you fail to reject the null hypothesis of the Breusch-Pagan test, then heteroscedasticity is not present and you can proceed to interpret the output of the original regression.
However, if you reject the null hypothesis, this means heteroscedasticity is present in the data. In this case, the standard errors that are shown in the output table of the regression may be unreliable.
There are a couple common ways that you can fix this issue, including:
1. Transform the response variable. You can try performing a transformation on the response variable. For example, you could use the log of the response variable instead of the original response variable. Typically taking the log of the response variable is an effective way of making heteroscedasticity go away. Another common transformation is to use the square root of the response variable.
2. Use weighted regression. This type of regression assigns a weight to each data point based on the variance of its fitted value. Essentially, this gives small weights to data points that have higher variances, which shrinks their squared residuals. When the proper weights are used, this can eliminate the problem of heteroscedasticity.
Hello
I tried to load the library for lmtest in R from the console. A comment came back to me there is no package called lmtest. On my computer the R-Stat package is loaded. I checked the table of contents under the R stats package and there is no lmtest. Could you help
Thank you!
Hello,
How do I apply the test when the model is a weighted linear regression? I want to check if the residuals are homoscedastic for a model using a weighting factor of 1/x^2 and another using 1/y^2. I tried using the unweighted model and assuming the weighting only in the bptest call, but it doesn’t seem to have worked. In addition, using Bartlett’s test, the residuals were homoscedastic. The residuals followed a normal distribution.
Regards
Hi Gustavo…Great question — a couple of subtle things are tripping you up.
First, if you fit the model *unweighted* and then try to “tell” `bptest` about the weights, you’ll get a misleading result. Heteroskedasticity tests use the *residuals from the fitted model*, so you want to test the residuals from the **weighted** fit, not from an ordinary least squares fit.
Here’s a clean way to do it in R.
1. Fit the weighted models and run BP/White-style tests on those fits
“`r
library(lmtest)
library(sandwich)
# weights ~ 1/x^2
fit_x <- lm(y ~ x + z, data = df, weights = I(1 / x^2)) bptest(fit_x, studentize = TRUE) # Breusch–Pagan bptest(fit_x, ~ fitted(fit_x) + I(fitted(fit_x)^2)) # White-type variant # If you still worry about heteroskedasticity, use robust SEs: coeftest(fit_x, vcov = vcovHC(fit_x, type = "HC3")) ``` 2. Check the *right* residuals For WLS, look at standardized or Pearson residuals, which account for the weights: ```r plot(fitted(fit_x), rstandard(fit_x)) # should look like constant spread plot(fitted(fit_x), resid(fit_x, type="pearson")) # Also try the scale–location plot: plot(fit_x, which = 3) ``` 3. About using 1/y^2 as weights Be careful here. Good weights come from prior knowledge of error variance, not from the response itself. Using 1/y^2 can induce bias because the weight depends on a noisy quantity on the left-hand side. If those weights reflect known measurement-error variances of y, fine; otherwise it’s safer to: * model the variance explicitly, or * estimate a variance–mean relationship from residuals and use iteratively reweighted least squares (IRLS). A principled alternative is `gls()` with a variance function: ```r library(nlme) # If Var(e_i) ∝ x^2, use varFixed(~ 1/x) fit_gls <- gls(y ~ x + z, data = df, weights = varFixed(~ 1/x)) # Compare to unweighted via likelihood ratio: anova(update(fit_gls, weights = NULL), fit_gls) ``` 4. Why Bartlett’s said “homoscedastic” Bartlett’s test is mainly for comparing group variances and is sensitive to normality assumptions. It’s not ideal when heteroskedasticity varies smoothly with a continuous predictor. Your residuals being normal is good news, but it doesn’t guarantee constant variance across x (or fitted values). 5. Quick checklist * Refit with the intended weights and test *that model’s* residuals. * Prefer BP/White tests and visual checks over Bartlett’s for continuous predictors. * Inspect standardized/pearson residuals vs fitted and vs the variable used in the weights (x). * If heteroskedasticity persists, keep the WLS fit but report robust (HC3) standard errors for inference. If you share a tiny reproducible snippet (formula, a couple of columns, and which extra predictors are in the model), I can tailor the variance model (e.g., whether to use \~x, \~log(x), or add a quadratic term in the BP test) and sanity-check whether 1/x^2 is doing what you expect.