0
$\begingroup$

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!

Edit: Added residuals plot Model 1 Model 2

Improve this question
$\endgroup$
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$
    – medium-dimensional
    Commented Aug 27, 2023 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
    Commented Aug 27, 2023 at 12:02

1 Answer 1

Reset to default
0
$\begingroup$

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
Improve this answer
$\endgroup$
2
  • $\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
    Commented Aug 27, 2023 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.
    Commented Aug 27, 2023 at 14:09

Your Answer

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

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