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

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)?

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.

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.