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

Estoy ejecutando una regresión OLS agrupada usando el paquete plm en R. Sin embargo, mi pregunta es más sobre estadísticas básicas, por lo que intento publicarlo aquí primero;)

Dado que mis resultados de regresión producen residuos heteroscedasticos, me gustaría intentar usar errores estándar robustos de heteroscedasticidad. Como resultado 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 analizar la importancia de las diferentes variables, me gustaría trazar la proporción de varianza explicada por cada variable independiente, por lo que necesito la suma de cuadrados correspondiente. Sin embargo, usando la función aov(), no sé cómo decirle a R que use errores estándar robustos.

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

Editar:

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

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

Editar:

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

El ANOVA en modelos de regresión lineal es equivalente a la prueba de Wald (y la prueba de razón de verosimilitud) de los modelos anidados correspondientes. Entonces, cuando desee realizar la prueba correspondiente utilizando errores estándar de coherencia a la heterocedasticidad (HC), esto no puede obtenerse a partir de una descomposición de las sumas de cuadrados, pero puede realizar la prueba de Wald utilizando una estimación de covarianza HC. Esta idea se utiliza tanto en Anova()y linearHypothesis()desde el carpaquete y coeftest()y waldtest()del lmtestpaquete. Los últimos tres también se pueden usar con plmobjetos.

Un ejemplo simple (aunque no muy interesante / significativo) es el siguiente. Utilizamos el modelo estándar de la ?plmpágina del manual y queremos llevar a cabo una prueba de Wald para determinar la importancia 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"))

Primero, veamos las pruebas marginales de Wald con errores estándar de 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 llevamos a cabo 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 ( mod0digamos) sin los dos coeficientes y luego llamar 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 prueba y el valor p calculados por linearHypothesis()y waldtest()son exactamente los mismos. Solo la interfaz y el formato de salida es algo diferente. En algunos casos, uno u otro es más simple o más intuitivo.

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

Esto se puede hacer con la Anovafunción en el carpaquete. Consulte esta excelente respuesta para obtener más detalles y una revisión de otras técnicas para tratar la heterocedasticidad en ANOVA.