В чем разница между этими двумя тестами Бреуша-язычества?

Используя R на некоторых данных и пытаясь определить , являются ли мои данные гетероскедастичными, я нашел две реализации теста Бреуша -Пагана: bptest (package lmtest) и ncvTest (package car). Однако они дают разные результаты. Какая разница между двумя? Когда вы должны использовать один или другой?

> 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 

Этот пример показывает, что согласно тестам мои данные в одном случае гетероскедастичны, а в другом - гомоскедастичны. Я нашел этот вопрос здесь, чтобы bptest мог быть изучен, а ncvTest - нет, однако, что это значит?

Ваше предположение верно, ncvTestвыполняет оригинальную версию теста Бреуша-Пагана. Это можно проверить, сравнив его с bptest(model, studentize = FALSE). (Как указывало @ Helix123, две функции также отличаются в других аспектах, таких как аргументы по умолчанию, следует проверить руководства по пакетам lmtestи получить carболее подробную информацию.)

$\xi^\ast$$\hat{\xi}$

Вот фрагмент кода, который демонстрирует то, что я только что написал (данные взяты из farawayпакета):

> 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

Что касается того, почему кто-то хочет изучать оригинальный тест АД, может помочь прямая цитата из статьи Р. Кенкера:

... Из этого анализа вытекают два вывода:

1. Асимптотическая сила теста Брейша и Пэгана чрезвычайно чувствительна к эксцессу распределения , и$\varepsilon$
2. Асимптотический размер теста верен только в частном случае гауссовского эксцесса.

Первый вывод расширен в работе Koenker and Bassett (1981), где предлагаются альтернативные, надежные тесты на гетероскедастичность. Последний вывод подразумевает, что уровни значимости, предложенные Брейшем и Паганом, будут правильными только в гауссовых условиях на . Поскольку такие условия обычно предполагаются на основе слепой веры и общеизвестно, что их трудно проверить, предлагается модификация теста Брейша и Пэгана, которая корректно «обучает» статистику теста и приводит к асимптотически корректным уровням значимости для достаточно большого класса распределений для .$\varepsilon$$\varepsilon$

Короче говоря, проверенный BP тест является более надежным, чем исходный.

В практическом плане ncvTestиспользует левую часть уравнения и bptestиспользует правую часть по умолчанию.

Это означает, что в случае Y ~ Xобоих тестов будут получены одинаковые результаты (относительно studentize = Fопции bptest). Но в многомерном анализе, например Y ~ X1 + X2, результаты будут отличаться. (Как указал @ Helix123)

Из справочного файла ncvTest : var.formula: «Односторонняя формула для дисперсии ошибок; если она не указана, дисперсия ошибок зависит от подгоночных значений ». Это означает, что по умолчанию будут использоваться подогнанные значения, но это также позволяет использовать линейную комбинацию независимых переменных (X1 + X2).

Из справочного файла bptest : varformula: «По умолчанию используются те же пояснительные переменные , что и в основной модели регрессии».

Продолжая тот же пример @Francis (данные stat500, из farawayпакета):

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

Тест BP, используя установленные значения:

> 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

Тест АД с использованием линейной комбинации предикторов:

> 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

«Опция линейной комбинации» позволяет исследовать гетероскедастичность, связанную с линейной зависимостью конкретной независимой переменной. Например, просто hwпеременная:

> 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

И наконец, как подытожил @Francis: «Короче говоря, тест BP с проверкой на прочность является более надежным, чем исходный», я обычно bptestиспользую studentize = TRUE(по умолчанию) и в varformula = ~ fitted.values(my.lm)качестве опций первоначальный подход к гомоскедастичности.