Comment obtenir une table ANOVA avec des erreurs standard robustes?

10

J'exécute une régression OLS groupée à l'aide du package plm dans R. Bien que ma question concerne davantage les statistiques de base, j'essaie donc de la publier ici en premier;)

Étant donné que mes résultats de régression produisent des résidus hétéroscédastiques, je voudrais essayer d'utiliser des erreurs standard robustes d'hétéroskédasticité. Comme résultat, coeftest(mod, vcov.=vcovHC(mod, type="HC0"))j'obtiens un tableau contenant des estimations, des erreurs types, des valeurs t et des valeurs p pour chaque variable indépendante, qui sont essentiellement mes résultats de régression "robustes".

Pour discuter de l'importance des différentes variables, je voudrais tracer la part de variance expliquée par chaque variable indépendante, j'ai donc besoin de la somme respective des carrés. Cependant, en utilisant la fonction aov(), je ne sais pas comment dire à R d'utiliser des erreurs standard robustes.

Maintenant, ma question est: comment puis-je obtenir la table / somme des carrés ANOVA qui fait référence à des erreurs standard robustes? Est-il possible de le calculer sur la base du tableau ANOVA à partir d'une régression avec des erreurs types normales?

Éditer:

En d'autres termes et sans tenir compte de mes problèmes R:

Si R 2 n'est pas affecté par l'utilisation d'erreurs types robustes, les contributions respectives à la variance expliquée par les différentes variables explicatives seront-elles également inchangées?2

Éditer:

Dans R, donne-t-on aov(mod)réellement une table ANOVA correcte pour un panelmodel (plm)?

Aki
la source

Réponses:

12

L'ANOVA dans les modèles de régression linéaire est équivalente au test de Wald (et au test du rapport de vraisemblance) des modèles imbriqués correspondants. Ainsi, lorsque vous souhaitez effectuer le test correspondant à l'aide d'erreurs standard à hétéroscédasticité cohérente (HC), cela ne peut pas être obtenu à partir d'une décomposition des sommes des carrés, mais vous pouvez effectuer le test de Wald à l'aide d'une estimation de covariance HC. Cette idée est utilisée à la fois Anova()et linearHypothesis()de l' caremballage et coeftest()et waldtest()du lmtestpaquet. Les trois derniers peuvent également être utilisés avec des plmobjets.

Un exemple simple (quoique peu intéressant / significatif) est le suivant. Nous utilisons le modèle standard de la ?plmpage de manuel et voulons effectuer un test de Wald pour la signification des deux log(pcap)et unemp. Nous avons besoin de ces packages:

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

Le modèle (sous l'alternative) est:

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

Tout d'abord, regardons les tests de Wald marginaux avec des erreurs standard HC pour tous les coefficients individuels:

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

Et puis nous effectuons un test de Wald pour les deux log(pcap)et 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

Alternativement, nous pouvons également ajuster le modèle sous l'hypothèse nulle ( mod0disons) sans les deux coefficients, puis appeler 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

La statistique de test et la valeur de p calculées par linearHypothesis()et waldtest()sont exactement les mêmes. Le formatage de l'interface et de la sortie est quelque peu différent. Dans certains cas, l'un ou l'autre est plus simple ou plus intuitif.

Remarque: Si vous fournissez une estimation de matrice de covariance (c.-à-d. Une matrice comme vocvHC(mod)) au lieu d'un estimateur de matrice de covariance (c.-à-d. Une fonction comme vocvHC), assurez-vous de fournir l'estimation de matrice de covariance HC du modèle sous l'alternative, c.-à-d. modèle non restreint.

Achim Zeileis
la source
1
Si je comprends bien, le test de Wald montre si l'inclusion de certaines variables est significative ou non. Mais existe-t-il une mesure de l' amélioration réelle du modèle, par exemple la somme des carrés?
Aki
Comment puis-je implémenter les erreurs standard HAC? J'ai essayé coeftest (mod, vcov = vcovHAC) mais j'ai reçu un message d'erreur disant "aucune méthode applicable pour 'estfun' appliquée à un objet de classe" c ('plm', 'panelmodel') ". Il semble que je doive calculer des poids ou une fonction d'estimation en premier. Quelle méthode recommanderiez-vous?
Aki
Bien que le plmpackage ait des méthodes pour le vcovHC()générique du sandwichpackage, il ne fournit pas de méthodes pour vcovHAC(). Au lieu de cela, il est plmlivré avec ses propres fonctions dédiées au calcul des matrices de covariance dans les modèles de panneaux qui peuvent également inclure une corrélation série. Voir vcovNW()ou vcovSCC()en emballage plm.
Achim Zeileis
Je vous remercie! Pour autant que je comprends, les deux fonctions se rapportent à un SE robuste à l'autocorrélation. Y a-t-il une fonction qui peut être utilisée pour les SE robustes à l'hétéroscédasticité et à l'autocorrélation?
Aki
Les trois fonctions ( vcovHAC, vcovNW, vcovSCC) sont _H_eteroskedasticity et _A_utocorrelation _C_onsistent ... Voilà ce que AHC défend.
Achim Zeileis
2

Cela peut être fait avec la Anovafonction dans le carpackage. Voir cette excellente réponse pour plus de détails et une revue d'autres techniques pour traiter l'hétéroskédasticité en ANOVA.

shadowtalker
la source
Je vous remercie. Le problème est que Anova () ne semble pas fonctionner avec des modèles de type plm (panelmodel).
Aki
@Aki si je ne me trompe pas, l'OLS groupé est équivalent à l'OLS ordinaire, basé sur ce qu'il dit dans la vignette: cran.r-project.org/web/packages/plm/vignettes/plm.pdf
shadowtalker
@Aki cependant, il semble que vous pourriez être intéressé par un modèle ANOVA plus riche. Il y a quelques exemples de R ici: stats.stackexchange.com/questions/3874/…
shadowtalker