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 car
pakietu i coeftest()
oraz waldtest()
z lmtest
pakietu. Tych ostatnich można również używać z plm
obiektami.
Prosty (choć niezbyt interesujący / znaczący) przykład jest następujący. Używamy standardowego modelu ze strony ?plm
podrę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ą ( mod0
powiedzmy) 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.
plm
pakiet zawiera metodyvcovHC()
ogólne zsandwich
pakietu, nie udostępnia metodvcovHAC()
. Zamiast tegoplm
dostarczany jest z własnymi dedykowanymi funkcjami do obliczania macierzy kowariancji w modelach panelowych, które potencjalnie obejmują również korelację szeregową. ZobaczvcovNW()
lubvcovSCC()
w paczceplm
.vcovHAC
,vcovNW
,vcovSCC
) są _H_eteroskedasticity i _A_utocorrelation _C_onsistent ... to co stoi na HAC.Można to zrobić za pomocą
Anova
funkcji wcar
pakiecie. Zobacz tę doskonałą odpowiedź, aby uzyskać więcej szczegółów i przegląd innych technik radzenia sobie z heteroskedastycznością w ANOVA.źródło