Jaka jest różnica między tymi dwoma testami Breuscha-Pagana?

Używając R na niektórych danych i próbując sprawdzić, czy moje dane są heteroscedastyczne, znalazłem dwie implementacje testu Breusch-Pagan, bptest (pakiet lmtest) i ncvTest (pakiet samochodów). Dają one jednak różne wyniki. Jaka jest różnica między nimi? Kiedy powinieneś wybrać jedno lub drugie?

> 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 

Ten przykład pokazuje, że zgodnie z testami moje dane są w jednym przypadku heteroscedastyczne, aw drugim przypadku homoscedastyczne. Znalazłem to pytanie tutaj, więc bptest może być uczony, a ncvTest może nie być, ale co to znaczy?

Twoje przypuszczenie jest prawidłowe, ncvTestwykonuje oryginalną wersję testu Breusch-Pagana. Można to faktycznie zweryfikować, porównując z tym bptest(model, studentize = FALSE). (Jak wskazał @ Helix123, dwie funkcje różnią się także innymi aspektami, takimi jak domyślne argumenty, należy sprawdzić instrukcje pakietów lmtesti caruzyskać więcej szczegółów.)

Studencki test Breuscha-Pagana został zaproponowany przez R. Koenkera w jego artykule z 1981 r. A Note on Studentizing a Test for Heteroscedasticity . Najbardziej oczywistą różnicą między nimi jest to, że używają różnych statystyk testowych. Mianowicie, niech$\xi^\ast$ będą statystycznymi badaniami testowymi i$\hat{\xi}$bądź oryginalny,

Oto fragment kodu, który pokazuje, co właśnie napisałem (dane pobrane z farawaypakietu):

> 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

Jeśli chodzi o to, dlaczego chcemy studiować oryginalny test BP, pomocny może być bezpośredni cytat z artykułu R. Koenkera:

... Z tej analizy wynikają dwa wnioski:

1. Asymptotyczna moc testu Breuscha i Pagana jest niezwykle wrażliwa na kurtozę dystrybucji , oraz$\varepsilon$
2. asymptotyczny rozmiar testu jest poprawny tylko w szczególnym przypadku kurtozy gaussowskiej.

Pierwszy wniosek został rozwinięty w Koenker i Bassett (1981), gdzie sugerowane są alternatywne, solidne testy heteroscedastyczności. Ten ostatni wniosek sugeruje, że poziomy istotności sugerowane przez Breuscha i Pagana będą poprawne tylko w warunkach gaussowskich na . Ponieważ takie warunki są ogólnie zakładane w ślepej wierze i są niezwykle trudne do zweryfikowania, sugerowana jest modyfikacja testu Breuscha i Pagana, która prawidłowo „uczą się” statystyki testu i prowadzi do asymptotycznie poprawnych poziomów istotności dla stosunkowo dużej klasy rozkładów dla .$\varepsilon$$\varepsilon$

Krótko mówiąc, studencki test BP jest bardziej niezawodny niż test oryginalny.

W praktyce ncvTestużywa lewej strony równania i bptestdomyślnie używa prawej.

Oznacza to, że w przypadku Y ~ X, oba testy dostarczą takie same wyniki (w odniesieniu do studentize = Fopcji bptest). Ale w analizie wielowymiarowej, takiej jak Y ~ X1 + X2, wyniki będą inne. (Jak wskazał @ Helix123)

Z pliku pomocy ncvTest : var.formula: „jednostronna formuła wariancji błędu; jeśli zostanie pominięta, wariancja błędu zależy od dopasowanych wartości ”. Co oznacza, że ​​domyślnie zostaną zastosowane dopasowane wartości, ale pozwala również na użycie liniowej kombinacji zmiennych niezależnych (X1 + X2).

Z pliku pomocy bptest : varformula: „Domyślnie pobierane są te same zmienne objaśniające , co w głównym modelu regresji”.

Kontynuując ten sam przykład @Francis (dane stat500, z farawaypakietu):

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

Test BP, przy użyciu dopasowanych wartości:

> 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

Test BP z wykorzystaniem liniowej kombinacji predyktorów:

> 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

„Opcja kombinacji liniowej” pozwala zbadać heteroskedastyczność związaną z liniową zależnością określonej zmiennej niezależnej. Na przykład tylko hwzmienna:

> 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

Wreszcie, jak podsumował @Francis: „Krótko mówiąc, studencki test BP jest bardziej niezawodny niż test oryginalny”, zwykle używam bptest, z studentize = TRUE(domyślnie) i varformula = ~ fitted.values(my.lm)jako opcjami, do wstępnego podejścia do homoskedastyczności.