What is the difference between these two Breusch-Pagan Tests?

Question

Using R on some data and trying to see whether or not my data is heteroscedastic, I've found two implementations of the Breusch-Pagan test, bptest (package lmtest) and ncvTest (package car). However, these produce different results. What is the difference between the two? When should you choose to use one or the other?

> model <- lm(y ~ x)
> bp <- bptest(model)
> bp
studentized Breusch-Pagan test

data:  model 
BP = 3.3596, df = 1, p-value = 0.06681

> ncvTest(model)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 3.858704    Df = 1     p = 0.04948855 

These example shows that according to the tests, my data is in one case heteroscedastic and in the other case homoscedastic. I did find this question here so bptest might be studentized and ncvTest might not be, however, what does this mean then?

Answer 1

Your guess is correct, ncvTest performs the original version of Breusch-Pagan test. This can actually be verified by comparing it to bptest(model, studentize = FALSE). (As @Helix123 pointed out, two functions also differ in other aspects such as default arguments, one should check package manuals of lmtest and car for more detail.)

The studentized Breusch-Pagan test was proposed by R. Koenker in his 1981 article A Note on Studentizing a Test for Heteroscedasticity. The most obvious difference of the two is that they use different test statistics. Namely, let $\xi^\ast$ be the studentized test statistics and $\hat{\xi}$ be the original one, $$\newcommand{\Var}{\operatorname{Var}}\hat{\xi}=\lambda\xi^\ast,\qquad\lambda=\frac{\Var(\varepsilon^2)}{2\Var(\varepsilon)^2}.$$

Here is a snippet of code that demonstrates what I just wrote (data taken from faraway package):

> mdl = lm(final ~ midterm, data = stat500)
> bptest(mdl)

    studentized Breusch-Pagan test

data:  mdl
BP = 0.86813, df = 1, p-value = 0.3515

> bptest(mdl, studentize = FALSE)

    Breusch-Pagan test

data:  mdl
BP = 0.67017, df = 1, p-value = 0.413

> ncvTest(mdl)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.6701721    Df = 1     p = 0.4129916 
> 
> n = nrow(stat500)
> e = residuals(mdl)
> bpmdl = lm(e^2 ~ midterm, data = stat500)
> lambda = (n - 1) / n * var(e^2) / (2 * ((n - 1) / n * var(e))^2)
> Studentized_bp = n * summary(bpmdl)$r.squared
> Original_bp = Studentized_bp * lambda
> 
> Studentized_bp
[1] 0.8681335
> Original_bp
[1] 0.6701721

As for why one wants to studentize the original BP test, a direct quote from R. Koenker's article may be helpful:

... Two conclusions emerge from this analysis:

1. The asymptotic power of the Breusch and Pagan test is extremely sensitive to the kurtosis of the distribution of $\varepsilon$, and
2. the asymptotic size of the test is correct only in special case of Gaussian kurtosis.

The former conclusion is expanded upon in Koenker and Bassett (1981) where alternative, robust tests for heteroscedasticity are suggested. The latter conclusion implies that the significance levels suggested by Breusch and Pagan will be correct only under Gaussian conditions on $\varepsilon$. Since such conditions are generally assumed on blind faith and are notoriously difficult to verify, a modification of the Breusch and Pagan test is suggested which correctly "studentise" the test statistic and leads to asymptotically correct significance levels for a reasonably large class of distributions for $\varepsilon$.

In short, the studentized BP test is more robust than the original one.

Answer 2

In practical terms, ncvTest uses the left-hand-side of the equation and bptest uses the right-hand-side, by default.

It means that in a case of Y ~ X, both tests will provide the same results (regarding the studentize = F option of bptest). But in a multivariate analysis such as Y ~ X1 + X2, the results will be different. (As @Helix123 pointed out)

From the help file of ncvTest : var.formula : "a one-sided formula for the error variance; if omitted, the error variance depends on the fitted values." Which means that, by default, the fitted values will be used, but it also allows to use a linear combination of the independent variables (X1 + X2).

From the help file of bptest : varformula : "By default the same explanatory variables are taken as in the main regression model."

Continuing the same example of @Francis (data stat500, from faraway package):

> mdl_t = lm(final ~ midterm + hw, data = stat500)

BP test, using fitted values:

> ncvTest(mdl_t) # Default

Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.6509135    Df = 1     p = 0.4197863 

> bptest(mdl_t, varformula = ~ fitted.values(mdl_t), studentize = F)

Breusch-Pagan test

data:  mdl_t
BP = 0.65091, df = 1, p-value = 0.4198

BP test, using a linear combination of predictors:

> ncvTest(mdl_t, var.formula = ~ midterm + hw)
Non-constant Variance Score Test 
Variance formula: ~ midterm + hw 
Chisquare = 0.7689743    Df = 2     p = 0.6807997 

> bptest(mdl_t, studentize = F) # Default

Breusch-Pagan test

data:  mdl_t
BP = 0.76897, df = 2, p-value = 0.6808

The "linear combination option" allows to investigate heteroskedasticity associated to linear dependence of a specific independent variable. For example, just the hw variable:

> ncvTest(mdl_t, var.formula = ~ hw)
Non-constant Variance Score Test 
Variance formula: ~ hw 
Chisquare = 0.04445789    Df = 1     p = 0.833004 

> bptest(mdl_t, varformula = ~ stat500$hw, studentize = F)

Breusch-Pagan test

data:  mdl_t
BP = 0.044458, df = 1, p-value = 0.833

Lastly, as @Francis summarized, "In short, the studentized BP test is more robust than the original one", I usually go with bptest, with studentize = TRUE (default) and varformula = ~ fitted.values(my.lm) as options, for an initial approach for homoskedasticity.