# Heteroskedasticity testing with Breusch-Pagan in R, two functions — very different answers?

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.

• Playing around a bit more with the two functions answered my question. car::ncvTest(), by default, fits residuals to fitted/estimated values from the regression model only. lmtest::bptest() fits the residuals to all terms of the regression model by default. Therefore, both models answer different questions. – Zage12 May 29 at 3:41