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.

  • $\begingroup$ 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. $\endgroup$ – Zage12 May 29 at 3:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.