Wie bekomme ich eine ANOVA-Tabelle mit robusten Standardfehlern?

10

Ich führe eine gepoolte OLS-Regression mit dem plm-Paket in R aus. Meine Frage bezieht sich jedoch eher auf grundlegende Statistiken. Deshalb versuche ich, sie zuerst hier zu veröffentlichen.

Da meine Regressionsergebnisse heteroskedastische Residuen ergeben, möchte ich versuchen, robuste Standardfehler der Heteroskedastizität zu verwenden. Als Ergebnis coeftest(mod, vcov.=vcovHC(mod, type="HC0"))erhalte ich eine Tabelle mit Schätzungen, Standardfehlern, t-Werten und p-Werten für jede unabhängige Variable, die im Grunde meine "robusten" Regressionsergebnisse sind.

Um die Wichtigkeit verschiedener Variablen zu diskutieren, möchte ich den Anteil der Varianz darstellen, der durch jede unabhängige Variable erklärt wird, daher benötige ich die jeweilige Summe der Quadrate. Mit der Funktion aov()weiß ich jedoch nicht, wie ich R anweisen soll, robuste Standardfehler zu verwenden.

Meine Frage lautet nun: Wie erhalte ich die ANOVA-Tabelle / Quadratsumme, die sich auf robuste Standardfehler bezieht? Ist es möglich, es basierend auf der ANOVA-Tabelle aus der Regression mit normalen Standardfehlern zu berechnen?

Bearbeiten:

Mit anderen Worten und ohne Rücksicht auf meine R-Probleme:

Wenn R nicht durch die Verwendung robuster Standardfehler beeinflusst wird, bleiben dann auch die jeweiligen Beiträge zur erklärten Varianz durch die verschiedenen erklärenden Variablen unverändert?2

Bearbeiten:

aov(mod)Gibt es in R tatsächlich eine korrekte ANOVA-Tabelle für ein Panelmodell (plm)?

Aki
quelle

Antworten:

12

Die ANOVA in linearen Regressionsmodellen entspricht dem Wald-Test (und dem Likelihood-Ratio-Test) der entsprechenden verschachtelten Modelle. Wenn Sie also den entsprechenden Test mit heteroskedastizitätskonsistenten (HC) Standardfehlern durchführen möchten, kann dies nicht aus einer Zerlegung der Quadratsummen erhalten werden, sondern Sie können den Wald-Test mit einer HC-Kovarianzschätzung durchführen. Diese Idee wird sowohl in Anova()als auch linearHypothesis()aus dem carPaket und coeftest()und waldtest()aus dem lmtestPaket verwendet. Die letzten drei können auch mit plmObjekten verwendet werden.

Ein einfaches (wenn auch nicht sehr interessantes / aussagekräftiges) Beispiel ist das folgende. Wir verwenden das Standardmodell aus der ?plmHandbuchseite und möchten einen Wald-Test für die Bedeutung von log(pcap)und durchführen unemp. Wir brauchen diese Pakete:

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

Das Modell (unter der Alternative) ist:

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

Betrachten wir zunächst die marginalen Wald-Tests mit HC-Standardfehlern für alle einzelnen Koeffizienten:

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

Und dann führen wir einen Wald-Test für beide durch log(pcap)und 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

Alternativ können wir das Modell auch unter die Nullhypothese ( mod0sagen wir) ohne die beiden Koeffizienten anpassen und dann aufrufen 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

Die von linearHypothesis()und berechnete Teststatistik und der p-Wert waldtest()sind genau gleich. Nur die Formatierung der Benutzeroberfläche und der Ausgabe ist etwas anders. In einigen Fällen ist das eine oder andere einfacher oder intuitiver.

Hinweis: Wenn Sie vocvHC(mod)anstelle eines Kovarianzmatrixschätzers (dh einer Funktion wie vocvHC) eine Kovarianzmatrixschätzung (dh eine Matrix wie ) angeben , stellen Sie sicher, dass Sie die HC-Kovarianzmatrixschätzung des Modells unter der Alternative, dh der nicht eingeschränktes Modell.

Achim Zeileis
quelle
1
Wenn ich es richtig verstehe, zeigt der Wald-Test, ob das Einbeziehen bestimmter Variablen signifikant ist oder nicht. Aber gibt es ein Maß dafür, um wie viel sie das Modell tatsächlich verbessern, wie z. B. die Summe der Quadrate?
Aki
Wie kann ich HAC-Standardfehler implementieren? Ich habe coeftest (mod, vcov = vcovHAC) versucht, aber eine Fehlermeldung erhalten, die besagt, dass "keine anwendbare Methode für 'estfun' auf ein Objekt der Klasse" c ('plm', 'panelmodel') angewendet wurde ". Es scheint, dass ich Gewichte berechnen muss oder eine Schätzfunktion zuerst. Welche Methode würden Sie empfehlen?
Aki
Während das plmPaket Methoden für das vcovHC()Generikum aus dem sandwichPaket enthält, bietet es keine Methoden für vcovHAC(). Wird stattdessen plmmit eigenen Funktionen zur Berechnung von Kovarianzmatrizen in Panel-Modellen geliefert, die möglicherweise auch eine serielle Korrelation enthalten. Siehe vcovNW()oder vcovSCC()im Paket plm.
Achim Zeileis
Vielen Dank! Soweit ich weiß, beziehen sich beide Funktionen auf eine autokorrelationssichere SE. Gibt es eine Funktion, die sowohl für heteroskedastische als auch für autokorrelationsresistente SE verwendet werden kann?
Aki
Alle drei Funktionen ( vcovHAC, vcovNW, vcovSCC) sind _H_eteroskedasticity und _A_utocorrelation _C_onsistent ... das ist , was HAC steht.
Achim Zeileis
2

Dies kann mit der AnovaFunktion im carPaket erfolgen. In dieser ausgezeichneten Antwort finden Sie weitere Einzelheiten und einen Überblick über andere Techniken zum Umgang mit Heteroskedastizität bei ANOVA.

Shadowtalker
quelle
Vielen Dank. Das Problem ist, dass Anova () nicht mit Modellen vom Typ plm (panelmodel) zu funktionieren scheint.
Aki
@Aki Wenn ich mich nicht irre, entspricht das gepoolte OLS dem einfachen OLS, basierend auf den Angaben
shadowtalker
@Aki jedoch klingt es so, als ob Sie an einem reichhaltigeren ANOVA-Modell interessiert sein könnten. Es gibt einige R Beispiele hier: stats.stackexchange.com/questions/3874/…
shadowtalker