Jak uzyskać tabelę ANOVA z niezawodnymi standardowymi błędami?

10

Korzystam z regresji puli OLS przy użyciu pakietu plm w R. Chociaż moje pytanie dotyczy bardziej podstawowych statystyk, więc najpierw postaram się je tutaj zamieścić;)

Ponieważ moje wyniki regresji dają reszty heteroskedastyczne, chciałbym spróbować użyć solidnych standardowych błędów heteroskedastycznych. W rezultacie coeftest(mod, vcov.=vcovHC(mod, type="HC0"))otrzymuję tabelę zawierającą oszacowania, błędy standardowe, wartości t i wartości p dla każdej zmiennej niezależnej, które w zasadzie są moimi „solidnymi” wynikami regresji.

Aby omówić znaczenie różnych zmiennych, chciałbym wykreślić udział wariancji wyjaśniony przez każdą zmienną niezależną, więc potrzebuję odpowiedniej sumy kwadratów. Jednak za pomocą funkcji aov()nie wiem, jak powiedzieć R, aby używał solidnych standardowych błędów.

Teraz moje pytanie brzmi: jak uzyskać tabelę ANOVA / sumę kwadratów, która odnosi się do solidnych błędów standardowych? Czy można to obliczyć na podstawie tabeli ANOVA z regresji z normalnymi błędami standardowymi?

Edytować:

Innymi słowy i ignorując moje problemy z R:

Gdy R 2 nie wpływa solidną błędów standardowych, również odpowiednie wkłady do objaśnione wariancji przez różne zmienne objaśniające niezmienione?2)

Edytować:

Czy w R aov(mod)rzeczywiście podaje prawidłową tabelę ANOVA dla modelu panelu (plm)?

Aki
źródło

Odpowiedzi:

12

ANOVA w modelach regresji liniowej jest równoważna testowi Walda (i testowi prawdopodobieństwa) odpowiednich modeli zagnieżdżonych. Więc jeśli chcesz przeprowadzić odpowiedni test z wykorzystaniem błędów standardowych zgodnych z heteroskedastycznością (HC), nie można tego uzyskać z rozkładu sum kwadratów, ale możesz przeprowadzić test Walda za pomocą oszacowania kowariancji HC. Pomysł ten jest wykorzystywany zarówno Anova()i linearHypothesis()od carpakietu i coeftest()oraz waldtest()z lmtestpakietu. Tych ostatnich można również używać z plmobiektami.

Prosty (choć niezbyt interesujący / znaczący) przykład jest następujący. Używamy standardowego modelu ze strony ?plmpodręcznika i chcemy przeprowadzić test Walda na znaczenie obu log(pcap)i unemp. Potrzebujemy tych pakietów:

library("plm")
library("sandwich")
library("car")
library("lmtest")

Model (w ramach alternatywy) to:

data("Produc", package = "plm")
mod <- plm(log(gsp) ~ log(pc) + log(emp) + log(pcap) + unemp,
  data = Produc, index = c("state", "year"))

Najpierw spójrzmy na marginalne testy Walda z błędami standardowymi HC dla wszystkich poszczególnych współczynników:

coeftest(mod, vcov = vcovHC)

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
log(pc)    0.2920069  0.0617425  4.7294 2.681e-06 ***
log(emp)   0.7681595  0.0816652  9.4062 < 2.2e-16 ***
log(pcap) -0.0261497  0.0603262 -0.4335   0.66480    
unemp     -0.0052977  0.0024958 -2.1226   0.03411 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Następnie przeprowadzamy test Walda dla obu log(pcap)i unemp:

linearHypothesis(mod, c("log(pcap)", "unemp"), vcov = vcovHC)

Linear hypothesis test

Hypothesis:
log(pcap) = 0
unemp = 0

Model 1: restricted model
Model 2: log(gsp) ~ log(pc) + log(emp) + log(pcap) + unemp

Note: Coefficient covariance matrix supplied.

  Res.Df Df  Chisq Pr(>Chisq)  
1    766                       
2    764  2 7.2934    0.02608 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Alternatywnie możemy również dopasować model pod hipotezą zerową ( mod0powiedzmy) bez dwóch współczynników, a następnie wywołać waldtest():

mod0 <- plm(log(gsp) ~ log(pc) + log(emp),
  data = Produc, index = c("state", "year"))
waldtest(mod0, mod, vcov = vcovHC)

Wald test

Model 1: log(gsp) ~ log(pc) + log(emp)
Model 2: log(gsp) ~ log(pc) + log(emp) + log(pcap) + unemp
  Res.Df Df  Chisq Pr(>Chisq)  
1    766                       
2    764  2 7.2934    0.02608 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Statystyka testowa i wartość p obliczone przez linearHypothesis()i waldtest()są dokładnie takie same. Po prostu interfejs i formatowanie wyjściowe jest nieco inne. W niektórych przypadkach jedno lub drugie jest prostsze lub bardziej intuicyjne.

Uwaga: Jeśli podasz oszacowanie macierzy kowariancji (tj. Macierz podobna vocvHC(mod)) zamiast estymatora macierzy kowariancji (tj. Funkcja podobna vocvHC), upewnij się, że podajesz oszacowanie macierzy kowariancji HC modelu zgodnie z alternatywą, tj. model nieograniczony.

Achim Zeileis
źródło
1
Jeśli dobrze to rozumiem, test Walda pokazuje, czy uwzględnienie niektórych zmiennych jest znaczące, czy nie. Ale czy istnieje miara tego, jak bardzo poprawiają one model, np. Suma kwadratów?
Aki
Jak mogę wdrożyć standardowe błędy HAC? Próbowałem coeftest (mod, vcov = vcovHAC), ale otrzymałem komunikat o błędzie „nie ma zastosowania do metody„ estfun ”zastosowanej do obiektu klasy„ c („plm”, „panelmodel”). Wydaje się, że muszę obliczyć wagi lub najpierw funkcja szacowania. Jaką metodę byś polecił?
Aki,
Chociaż plmpakiet zawiera metody vcovHC()ogólne z sandwichpakietu, nie udostępnia metod vcovHAC(). Zamiast tego plmdostarczany jest z własnymi dedykowanymi funkcjami do obliczania macierzy kowariancji w modelach panelowych, które potencjalnie obejmują również korelację szeregową. Zobacz vcovNW()lub vcovSCC()w paczce plm.
Achim Zeileis,
Dziękuję Ci! O ile rozumiem, obie funkcje odnoszą się do SE odpornego na autokorelację. Czy jest jakaś funkcja, która może być użyta do SE odpornego na heteroscedastyczność i autokorelację?
Aki
Wszystkie trzy funkcje ( vcovHAC, vcovNW, vcovSCC) są _H_eteroskedasticity i _A_utocorrelation _C_onsistent ... to co stoi na HAC.
Achim Zeileis
2

Można to zrobić za pomocą Anovafunkcji w carpakiecie. Zobacz tę doskonałą odpowiedź, aby uzyskać więcej szczegółów i przegląd innych technik radzenia sobie z heteroskedastycznością w ANOVA.

Shadowtalker
źródło
Dziękuję Ci. Problem polega na tym, że Anova () nie działa z modelami typu plm (panelmodel).
Aki
@ Aki, jeśli się nie mylę, połączony OLS jest równoważny zwykłemu OLS, w oparciu o to, co napisano w winiecie: cran.r-project.org/web/packages/plm/vignettes/plm.pdf
shadowtalker
@Aki jednak wydaje się, że możesz być zainteresowany bogatszym modelem ANOVA. Oto kilka przykładów R: stats.stackexchange.com/questions/3874/...
shadowtalker