¿Cómo obtener una tabla de ANOVA con errores estándar robustos?

Question

Estoy ejecutando una regresión OLS agrupada usando el paquete plm en R. Aunque, mi pregunta es más sobre estadística básica, así que intento publicarla aquí primero ;)

Dado que los resultados de mi regresión arrojan residuos heterocedásticos, me gustaría intentar utilizar errores estándar robustos a la heterocedasticidad. Como resultado de coeftest(mod, vcov.=vcovHC(mod, type="HC0")) Obtengo una tabla que contiene estimaciones, errores estándar, valores t y valores p para cada variable independiente, que básicamente son mis resultados de regresión "robustos".

Para discutir la importancia de las diferentes variables me gustaría trazar la parte de la varianza explicada por cada variable independiente, por lo que necesito la respectiva suma de cuadrados. Sin embargo, utilizando la función aov() No sé cómo decirle a R que utilice errores estándar robustos.

Ahora mi pregunta es: ¿Cómo obtengo la tabla ANOVA/suma de cuadrados que se refiere a los errores estándar robustos? ¿Es posible calcularla a partir de la tabla de ANOVA de la regresión con errores estándar normales?

Editar:

En otras palabras y sin tener en cuenta mis problemas con la R:

Si R $^2$ no se ve afectado por el uso de errores estándar robustos, ¿también se mantendrán las respectivas contribuciones a la varianza explicada por las diferentes variables explicativas?

Editar:

En R, ¿se aov(mod) ¿realmente da una tabla ANOVA correcta para un modelo de panel (plm)?

Answer 1

El ANOVA en los modelos de regresión lineal es equivalente a la prueba de Wald (y a la prueba de razón de verosimilitud) de los correspondientes modelos anidados. Por lo tanto, cuando se quiere realizar la prueba correspondiente utilizando errores estándar consistentes en la heteroscedasticidad (HC), ésta no se puede obtener a partir de una descomposición de las sumas de los cuadrados, pero se puede llevar a cabo la prueba de Wald utilizando una estimación de la covarianza HC. Esta idea se utiliza tanto en Anova() y linearHypothesis() de la car paquete y coeftest() y waldtest() de la lmtest paquete. Los tres últimos también pueden utilizarse con plm objetos.

Un ejemplo sencillo (aunque no muy interesante o significativo) es el siguiente. Utilizamos el modelo estándar del ?plm página del manual y quiere realizar una prueba de Wald para la significación de ambos log(pcap) y unemp . Necesitamos estos paquetes:

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

El modelo (bajo la alternativa) es:

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

En primer lugar, veamos las pruebas de Wald marginales con errores estándar HC para todos los coeficientes individuales:

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

Y luego realizamos una prueba de Wald para ambos log(pcap) y 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

Alternativamente, también podemos ajustar el modelo bajo la hipótesis nula ( mod0 digamos) sin los dos coeficientes y luego llamar 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

La estadística de la prueba y el valor p calculados por linearHypothesis() y waldtest() es exactamente lo mismo. Sólo la interfaz y el formato de salida son algo diferentes. En algunos casos uno u otro es más sencillo o intuitivo.

Nota: Si proporciona una estimación de la matriz de covarianza (es decir, una matriz como vocvHC(mod) ) en lugar de un estimador de la matriz de covarianza (es decir, una función como vocvHC ), asegúrese de suministrar la estimación de la matriz de covarianza HC del modelo bajo la alternativa, es decir, el modelo no restringido.

Answer 2

Esto se puede hacer con el Anova en la función car paquete. Véase esta excelente respuesta para más detalles y una revisión de otras técnicas para tratar la heteroscedasticidad en el ANOVA.