Agrupación de errores estándar en R (manualmente o en plm)

Question

Estoy tratando de entender el error estándar "clustering" y cómo ejecutarlo en R (es trivial en Stata). En R no he tenido éxito utilizando plm o escribir mi propia función. Utilizaré el diamonds datos del ggplot2 paquete.

Puedo hacer efectos fijos con variables ficticias

> library(plyr)
> library(ggplot2)
> library(lmtest)
> library(sandwich)
> # with dummies to create fixed effects
> fe.lsdv <- lm(price ~ carat + factor(cut) + 0, data = diamonds)
> ct.lsdv <- coeftest(fe.lsdv, vcov. = vcovHC)
> ct.lsdv

t test of coefficients:

                      Estimate Std. Error  t value  Pr(>|t|)    
carat                 7871.082     24.892  316.207 < 2.2e-16 ***
factor(cut)Fair      -3875.470     51.190  -75.707 < 2.2e-16 ***
factor(cut)Good      -2755.138     26.570 -103.692 < 2.2e-16 ***
factor(cut)Very Good -2365.334     20.548 -115.111 < 2.2e-16 ***
factor(cut)Premium   -2436.393     21.172 -115.075 < 2.2e-16 ***
factor(cut)Ideal     -2074.546     16.092 -128.920 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

o bien desindexando los lados izquierdo y derecho (aquí no hay regresores invariantes en el tiempo) y corrigiendo los grados de libertad.

> # by demeaning with degrees of freedom correction
> diamonds <- ddply(diamonds, .(cut), transform, price.dm = price - mean(price), carat.dm = carat  .... [TRUNCATED] 
> fe.dm <- lm(price.dm ~ carat.dm + 0, data = diamonds)
> ct.dm <- coeftest(fe.dm, vcov. = vcovHC, df = nrow(diamonds) - 1 - 5)
> ct.dm

t test of coefficients:

         Estimate Std. Error t value  Pr(>|t|)    
carat.dm 7871.082     24.888  316.26 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

No puedo replicar estos resultados con plm porque no tengo un índice de "tiempo" (es decir, no se trata realmente de un panel, sino de conglomerados que podrían tener un sesgo común en sus términos de error).

> plm.temp <- plm(price ~ carat, data = diamonds, index = "cut")
duplicate couples (time-id)
Error in pdim.default(index[[1]], index[[2]]) : 

También intenté codificar mi propia matriz de covarianza con error estándar agrupado utilizando la explicación de Stata sobre su cluster opción ( explicada aquí ), que consiste en resolver $$\hat V_{cluster} = (X'X)^{-1} \left( \sum_{j=1}^{n_c} u_j'u_j \right) (X'X)^{-1}$$ donde $u_j = \sum_{cluster~j} e_i * x_i$ , $n_c$ si el número de clusters, $e_i$ es el residuo del $i^{th}$ observación y $x_i$ es el vector de filas de los predictores, incluida la constante (esto también aparece como ecuación (7.22) en la obra de Wooldridge Datos transversales y de panel ). Pero el siguiente código da matrices de covarianza muy grandes. ¿Son estos valores muy grandes dado el pequeño número de clusters que tengo? Dado que no puedo obtener plm para hacer clusters en un factor, no estoy seguro de cómo evaluar mi código.

> # with cluster robust se
> lm.temp <- lm(price ~ carat + factor(cut) + 0, data = diamonds)
> 
> # using the model that Stata uses
> stata.clustering <- function(x, clu, res) {
+     x <- as.matrix(x)
+     clu <- as.vector(clu)
+     res <- as.vector(res)
+     fac <- unique(clu)
+     num.fac <- length(fac)
+     num.reg <- ncol(x)
+     u <- matrix(NA, nrow = num.fac, ncol = num.reg)
+     meat <- matrix(NA, nrow = num.reg, ncol = num.reg)
+     
+     # outer terms (X'X)^-1
+     outer <- solve(t(x) %*% x)
+ 
+     # inner term sum_j u_j'u_j where u_j = sum_i e_i * x_i
+     for (i in seq(num.fac)) {
+         index.loop <- clu == fac[i]
+         res.loop <- res[index.loop]
+         x.loop <- x[clu == fac[i], ]
+         u[i, ] <- as.vector(colSums(res.loop * x.loop))
+     }
+     inner <- t(u) %*% u
+ 
+     # 
+     V <- outer %*% inner %*% outer
+     return(V)
+ }
> x.temp <- data.frame(const = 1, diamonds[, "carat"])
> summary(lm.temp)

Call:
lm(formula = price ~ carat + factor(cut) + 0, data = diamonds)

Residuals:
     Min       1Q   Median       3Q      Max 
-17540.7   -791.6    -37.6    522.1  12721.4 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
carat                 7871.08      13.98   563.0   <2e-16 ***
factor(cut)Fair      -3875.47      40.41   -95.9   <2e-16 ***
factor(cut)Good      -2755.14      24.63  -111.9   <2e-16 ***
factor(cut)Very Good -2365.33      17.78  -133.0   <2e-16 ***
factor(cut)Premium   -2436.39      17.92  -136.0   <2e-16 ***
factor(cut)Ideal     -2074.55      14.23  -145.8   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1511 on 53934 degrees of freedom
Multiple R-squared: 0.9272, Adjusted R-squared: 0.9272 
F-statistic: 1.145e+05 on 6 and 53934 DF,  p-value: < 2.2e-16 

> stata.clustering(x = x.temp, clu = diamonds$cut, res = lm.temp$residuals)
                        const diamonds....carat..
const                11352.64           -14227.44
diamonds....carat.. -14227.44            17830.22

¿Se puede hacer esto en R? Es una técnica bastante común en econometría (hay un breve tutorial en esta conferencia ), pero no puedo resolverlo en R. ¡Gracias!

Answer 1

Para los errores estándar de los blancos agrupados por grupo con el plm marco de trabajo intentar

coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0",cluster="group"))

donde model.plm es un modelo plm.

Véase también este enlace

http://www.inside-r.org/packages/cran/plm/docs/vcovHC o la documentación del paquete plm

EDITAR:

Para la agrupación bidireccional (por ejemplo, grupo y tiempo), véase el siguiente enlace:

http://people.su.se/~ma/clustering.pdf

La agrupación y otra información, especialmente para Stata, se puede encontrar aquí:

http://www.kellogg.northwestern.edu/faculty/petersen/htm/papers/se/se_programming.htm

Answer 2

Creo que se produce cierta confusión al llamar "números imaginarios" a los números que sólo tienen parte imaginaria. Creo que si los matemáticos tuvieran otro nombre para darles (pero realmente no estoy pensando en un nombre mejor), habría menos confusión.

ACTUALIZACIÓN: Incluso una persona tan brillante como Gauss piensa que "imaginario" es una denominación terrible:

Los números imaginarios son reales [Parte 1: Introducción]

Answer 3

La forma más sencilla de calcular los errores estándar agrupados en R es utilizar la función de resumen modificada.

lm.object <- lm(y ~ x, data = data)
summary(lm.object, cluster=c("c"))

Hay un excelente Correo electrónico: sobre la agrupación en el marco de lm. El sitio también proporciona la función de resumen modificada para la agrupación en uno y dos sentidos. Puede encontrar la función y el tutorial aquí .

Answer 4

Si no tiene un time índice, no necesitas uno: plm añadirá uno ficticio por sí mismo, y no se utilizará a menos que lo pidas. Así que esta llamada debería funcionar :

> x <- plm(price ~ carat, data = diamonds, index = "cut")
 Error in pdim.default(index[[1]], index[[2]]) : 
  duplicate couples (time-id) 

Excepto que no lo hace, lo que sugiere que has dado con un error en plm . (Este error ya ha sido corregido en el SVN. Puede instalar la versión de desarrollo desde aquí .)

Pero como se trataría de una ficción time índice de todos modos, podemos crearlo nosotros mismos:

diamonds$ftime <- 1:NROW(diamonds)  ##fake time

Ahora esto funciona:

x <- plm(price ~ carat, data = diamonds, index = c("cut", "ftime"))
coeftest(x, vcov.=vcovHC)
## 
## t test of coefficients:
## 
##       Estimate Std. Error t value  Pr(>|t|)    
## carat  7871.08     138.44  56.856 < 2.2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Nota importante : vcovHC.plm() en plm estimará por defecto los SEs de Arellano agrupados por grupos. Lo cual es diferentes de lo que vcovHC.lm() en sandwich estimará (por ejemplo, el vcovHC SEs en la pregunta original), es decir, SEs consistentes con la heteroscedasticidad sin agrupación.

---

Un enfoque distinto es atenerse a lm regresiones de variables ficticias y el multiwayvcov paquete.

library("multiwayvcov")
fe.lsdv <- lm(price ~ carat + factor(cut) + 0, data = diamonds)
coeftest(fe.lsdv, vcov.= function(y) cluster.vcov(y, ~ cut, df_correction = FALSE))
## 
## t test of coefficients:
## 
##                      Estimate Std. Error t value  Pr(>|t|)    
## carat                 7871.08     138.44  56.856 < 2.2e-16 ***
## factor(cut)Fair      -3875.47     144.83 -26.759 < 2.2e-16 ***
## factor(cut)Good      -2755.14     117.56 -23.436 < 2.2e-16 ***
## factor(cut)Very Good -2365.33     111.63 -21.188 < 2.2e-16 ***
## factor(cut)Premium   -2436.39     123.48 -19.731 < 2.2e-16 ***
## factor(cut)Ideal     -2074.55      97.30 -21.321 < 2.2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

En ambos casos se obtendrán los SEs de Arellano (1987) con agrupación por grupos. El multiwayvcov es una evolución directa y significativa de las funciones de agrupación originales de Arai.

También se puede observar la matriz de varianza-covarianza resultante de ambos enfoques, que arroja la misma estimación de varianza para carat :

vcov.plm <- vcovHC(x)
vcov.lsdv <- cluster.vcov(fe.lsdv, ~ cut, df_correction = FALSE)
vcov.plm
##          carat
## carat 19165.28
diag(vcov.lsdv)
##                carat      factor(cut)Fair      factor(cut)Good factor(cut)Very Good   factor(cut)Premium     factor(cut)Ideal 
##            19165.283            20974.522            13820.365            12462.243            15247.584             9467.263