I am currently testing two different functions in R to determine heteroskedasticity in a regression model with 4 predictors, 2 continuous and 2 categorical. One function is telling me there is a high chance that the data is heteroskedastic, while the other is telling me it is most likely not. Why exactly are these two functions giving very different results?
The results of bptest() and ncvTest are as follows:
> bptest(model) studentized Breusch-Pagan test data: model BP = 18.613, df = 4, p-value = 0.0009363 > bptest(model, studentize=F) # For illustrative purposes Breusch-Pagan test data: model BP = 19.574, df = 4, p-value = 0.000606 > ncvTest(model) Non-constant Variance Score Test Variance formula: ~ fitted.values Chisquare = 0.001700664, Df = 1, p = 0.96711
As noted here and here, I am aware that the former defaults to the "studentized" version of the Breusch-Pagan test, while the former does not. However, in the example given by Francis here, the values between bptest(studentize=F) and ncvTest are identical, however in my case they wildly differ.
When I performed the manual calculation as suggested in Francis's answer to the previous question, the values I got were similar to those performed using the bptest function.