# Order of independent and dependent variable in Breusch-Pagan test

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)
bptest(model1)
model_2 <- lm(Dob~GERD_HRQL, data=myresearch)
bptest(model2)


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!

• 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. Commented Aug 27, 2023 at 10:54
• @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. Commented Aug 27, 2023 at 12:02

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,

library(lmtest)

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

lmtest::bptest(fit1)$$statistic BP 2.916375 lmtest::bptest(fit2)$$statistic
BP
0.0455314

# 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

• 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. Commented Aug 27, 2023 at 13:06
• 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. Commented Aug 27, 2023 at 14:09