Como obter a tabela ANOVA com erros padrão robustos?

10

Estou executando uma regressão OLS em pool usando o pacote plm em R. No entanto, minha pergunta é mais sobre estatísticas básicas, então tento publicá-la aqui primeiro;)

Como meus resultados de regressão produzem resíduos heterocedásticos, eu gostaria de tentar usar erros padrão robustos de heterocedasticidade. Como resultado coeftest(mod, vcov.=vcovHC(mod, type="HC0")), recebo uma tabela contendo estimativas, erros padrão, valores t e valores p para cada variável independente, que basicamente são meus resultados de regressão "robustos".

Para discutir a importância de diferentes variáveis, gostaria de plotar a parcela de variação explicada por cada variável independente, portanto, preciso da soma dos quadrados respectiva. No entanto, usando a função aov(), não sei como dizer ao R para usar erros padrão robustos.

Agora, minha pergunta é: Como obtenho a tabela ANOVA / soma dos quadrados que se refere a erros padrão robustos? É possível calculá-lo com base na tabela ANOVA a partir de regressão com erros padrão normais?

Editar:

Em outras palavras, e desconsiderando meus problemas de R:

Se R não for afetado pelo uso de erros padrão robustos, as respectivas contribuições para a variância explicada pelas diferentes variáveis ​​explicativas serão inalteradas?2

Editar:

Em R, aov(mod)realmente fornece uma tabela ANOVA correta para um modelo de painel (plm)?

Aki
fonte

Respostas:

12

A ANOVA nos modelos de regressão linear é equivalente ao teste de Wald (e o teste da razão de verossimilhança) dos modelos aninhados correspondentes. Portanto, quando você deseja realizar o teste correspondente usando erros padrão consistentes em heterocedasticidade (HC), isso não pode ser obtido a partir de uma decomposição das somas de quadrados, mas você pode executar o teste de Wald usando uma estimativa de covariância de HC. Esta ideia é usado em ambos Anova()e linearHypothesis()do carpacote e coeftest()e waldtest()do lmtestpacote. Os três últimos também podem ser usados ​​com plmobjetos.

Um exemplo simples (embora não muito interessante / significativo) é o seguinte. Usamos o modelo padrão da ?plmpágina de manual e queremos realizar um teste de Wald para a significância de ambos log(pcap)e unemp. Precisamos destes pacotes:

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

O modelo (sob a alternativa) é:

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

Primeiro, vejamos os testes marginais de Wald com erros padrão de HC para todos os coeficientes individuais:

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

E então realizamos um teste de Wald para ambos log(pcap)e 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

Como alternativa, também podemos ajustar o modelo sob a hipótese nula ( mod0digamos) sem os dois coeficientes e depois chamar 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

A estatística do teste e o valor p calculados linearHypothesis()e waldtest()são exatamente os mesmos. Apenas a interface e a formatação da saída são um pouco diferentes. Em alguns casos, um ou outro é mais simples ou mais intuitivo.

Nota: Se você fornecer uma estimativa da matriz de covariância (ou seja, semelhante a uma matriz vocvHC(mod)) em vez de um estimador de matriz de covariância (ou seja, uma função semelhante vocvHC), certifique-se de fornecer a estimativa da matriz de covariância HC do modelo sob a alternativa, ou seja, modelo não restrito.

Achim Zeileis
fonte
11
Se eu entendi direito, o teste de Wald mostra se a inclusão de determinadas variáveis ​​é significativa ou não. Mas há uma medida de quanto eles realmente melhoram o modelo, como por exemplo, soma dos quadrados?
Aki
Como posso implementar erros padrão do HAC? Eu tentei coeftest (mod, vcov = vcovHAC), mas recebi uma mensagem de erro dizendo "nenhum método aplicável para 'estfun' aplicado a um objeto da classe" c ('plm', 'panelmodel') ". Parece que preciso calcular pesos . ou uma função de estimativa primeira Qual o método que você recomendaria?
Aki
Embora o plmpacote possua métodos para o vcovHC()genérico do sandwichpacote, ele não fornece métodos para vcovHAC(). Em vez disso, é plmfornecido com suas próprias funções dedicadas para calcular matrizes de covariância em modelos de painel que potencialmente também incluem correlação serial. Veja vcovNW()ou vcovSCC()no pacote plm.
Achim Zeileis
Obrigado! Tanto quanto eu entendo, ambas as funções estão relacionadas ao SE robusto de autocorrelação. Existe alguma função que possa ser usada para o SE robusto de heterocedasticidade e autocorrelação?
Aki
Todas as três funções ( vcovHAC, vcovNW, vcovSCC) são _H_eteroskedasticity e _A_utocorrelation _C_onsistent ... Isso é o que HAC defende.
Achim Zeileis
2

Isso pode ser feito com a Anovafunção no carpacote. Veja esta excelente resposta para obter mais detalhes e uma revisão de outras técnicas para lidar com a heterocedasticidade na ANOVA.

shadowtalker
fonte
Obrigado. O problema é que o Anova () parece não funcionar com modelos do tipo plm (panelmodel).
Aki
@Aki, se não me engano, o OLS agrupado é equivalente ao OLS comum, com base no que diz na vinheta: cran.r-project.org/web/packages/plm/vignettes/plm.pdf
shadowtalker
@Aki, no entanto, parece que você pode estar interessado em um modelo ANOVA mais rico. Existem alguns exemplos de R aqui: stats.stackexchange.com/questions/3874/…
shadowtalker 6/15