I've been working on my master thesis, and I've been using R to test heteroskedasticity with Breusch-Pagan test. The code is rather simple:

model_1 <- lm(GERD_HRQL~Dob, data=myresearch)
model_2 <- lm(Dob~GERD_HRQL, data=myresearch)

And the results are, respectively:

model_1: BP = 0.52928, df = 1, p-value = 0.4669
model_2: BP = 4.6722, df = 1, p-value = 0.03065

And this really surprised me, because I thought that for a Pearson correlation with one dependent variable and one independent variable, the order in which one defines the variables should not affect the results of the Breusch-Pagan test for heteroscedasticity.

Could somebody please explain why those results differ? I don't know which value to report. If there is heteroscedasticity when I calculate a Pearson correlation, I would use a wild bootstrap method.

Thank you!

Edit: Added residuals plot Model 1 Model 2

  • $\begingroup$ Welcome to CV! Not directly related to your question, but it might be interesting to see how residual plots look like for the fitted models. May be you could add the residual plots (plot(model_x)) in the question (if possible)? It might also be helpful to read the description of ?bptest in R. $\endgroup$ Aug 27 at 10:54
  • 1
    $\begingroup$ @medium-dimensional Thank you. I have added a residual plots. I've the description in R for BP test, however unfortunately it didn't help. $\endgroup$
    – daniele
    Aug 27 at 12:02

1 Answer 1


As @medium-dimensional pointed out, the description of bptest tells us that "the Breusch-Pagan test fits a linear regression model to the residuals of a linear regression model..." (see also wiki page). However, residuals (hence, test statistic) will be different if you switch outcome and predictor. For example,


fit1 <- lm(mpg ~ disp, mtcars)
fit2 <- lm(disp ~ mpg, mtcars)


# Using auxiliary regression models (e.g., "regressing the squared residuals")
aux_fit1 <- lm(residuals(fit1)^2 ~ disp, mtcars)
aux_fit2 <- lm(residuals(fit2)^2 ~ mpg, mtcars)

nrow(mtcars) * summary(aux_fit1)$r.squared
[1] 2.916375
nrow(mtcars) * summary(aux_fit2)$r.squared  
[1] 0.0455314
  • $\begingroup$ Thank you. The problem is, I want to calculate intercorrelations between variables that I will in later analysis use as dependent variables. In regards to these variables, for the purpose of my thesis it's not important which variable predicts the other – the only thing that’s Important is the strength of their association. $\endgroup$
    – daniele
    Aug 27 at 13:06
  • 2
    $\begingroup$ No problem, I understand your issue but that's not what you asked. You asked why the bptest results differ when you switched the outcome and the predictor, and I tried to answer. So I guess, in your case, it does not make much sense to use bptest at all. We calculate residuals by taking the difference between observed and predicted values but you say it does not matter which variable predicts the other. So, maybe you have a different question in mind. In any case, before going further, I recommend you to check Nick Cox's answer here. $\endgroup$
    – T.E.G.
    Aug 27 at 14:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.