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

33

Estoy tratando de entender el error estándar "agrupamiento" y cómo ejecutarlo en R (es trivial en Stata). En RI no he tenido éxito usando plmo escribiendo mi propia función. Usaré los diamondsdatos del ggplot2paquete.

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 desanimando los lados izquierdo y derecho (no hay regresores invariables en el tiempo aquí) 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 plm, porque no tengo un índice de "tiempo" (es decir, esto no es realmente un panel, solo clústeres 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 de su clusteropción ( explicada aquí ), que es resolver donde , si el número de grupos, es el residuo para la observación y es el vector fila de predictores, incluida la constante (esto también aparece como la ecuación (7.22) en la sección transversal de Wooldridge y los datos del panel

{\hat{V}}_{do l tu s t mi r} = (X^{'} X)^{- 1} (\sum_{j = 1}^{{norte}_{do}} {tu}_{j}^{'} {tu}_{j}) (X^{'} X)^{- 1}

$\hat V_{cluster} = (X'X)^{-1} \left( \sum_{j=1}^{n_c} u_j'u_j \right) (X'X)^{-1}$

u_{j} = \sum_{c l u s t e r j} e_{i} * x_{i}

$u_j = \sum_{cluster~j} e_i * x_i$

n_{c}

$n_c$

e_{i}

$e_i$

i^{t h}

$i^{th}$

x_{i}

$x_i$ ) Pero el siguiente código proporciona matrices de covarianza muy grandes. ¿Son estos valores muy grandes dada la pequeña cantidad de grupos que tengo? Dado que no puedo plmhacer clusters en un factor, no estoy seguro de cómo comparar 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 entenderlo en R. ¡Gracias!

r panel-data standard-error fixed-effects-model clustered-standard-errors Richard Herron
fuente

77

@ricardh, maldiga a todos los economistas por no verificar si el término que quieren usar ya está en las estadísticas. Clúster en este contexto significa grupo y no tiene ninguna relación con el análisis de conglomerados, por eso rseek le dio resultados no relacionados. Por lo tanto, eliminé la etiqueta de agrupación. Para el análisis de datos del panel, consulte el paquete plm . Tiene una bonita viñeta, por lo que puede encontrar lo que desea. En cuanto a su pregunta, no está claro lo que quiere. ¿Dentro de los errores estándar del grupo?

mpiktas

@ricardh, ayudaría mucho si pudiera vincular a algún manual de Stata donde clusterse explica esta opción. Estoy seguro de que sería posible replicar en R.

mpiktas

2

+1 por ese comentario. Los economistas colonizan la terminología como locos. Aunque a veces es difícil elegir al villano. Me tomó un tiempo, por ejemplo, hasta que me di cuenta de que factorno tenía nada que ver factanalsino que se refería a variables categorizadas. Sin embargo, el clúster en R se refiere al análisis de clúster, k-means es simplemente EL método de partición: statmethods.net/advstats/cluster.html . No entiendo tu pregunta, pero también supongo que el clúster no tiene nada que ver con eso.

hans0l0

@mpiktas, @ ran2 - ¡Gracias! Espero haber aclarado la pregunta. En resumen, ¿por qué existe la "agrupación de errores estándar" si solo se trata de efectos fijos, que ya existían?

Richard Herron

1

La función cluster.vcov en el paquete "multiwayvcov" hace lo que está buscando.

29

Para los errores estándar blancos agrupados por grupo con el plmmarco, intente

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

donde model.plmes un modelo plm

Vea 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), consulte el siguiente enlace:

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

Aquí hay otra guía útil para el paquete plm específicamente que explica las diferentes opciones para los errores estándar agrupados:

http://www.princeton.edu/~otorres/Panel101R.pdf

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

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

EDITAR 2:

Aquí hay ejemplos que comparan R y stata: http://www.richard-bluhm.com/clustered-ses-in-r-and-stata-2/

Además, multiwayvcovpuede ser útil. Esta publicación proporciona una descripción general útil: http://rforpublichealth.blogspot.dk/2014/10/easy-clustered-standard-errors-in-r.html

De la documentación:

library(multiwayvcov)
library(lmtest)
data(petersen)
m1 <- lm(y ~ x, data = petersen)

# Cluster by firm
vcov_firm <- cluster.vcov(m1, petersen$firmid)
coeftest(m1, vcov_firm)
# Cluster by year
vcov_year <- cluster.vcov(m1, petersen$year)
coeftest(m1, vcov_year)
# Cluster by year using a formula
vcov_year_formula <- cluster.vcov(m1, ~ year)
coeftest(m1, vcov_year_formula)

# Double cluster by firm and year
vcov_both <- cluster.vcov(m1, cbind(petersen$firmid, petersen$year))
coeftest(m1, vcov_both)
# Double cluster by firm and year using a formula
vcov_both_formula <- cluster.vcov(m1, ~ firmid + year)
coeftest(m1, vcov_both_formula)

velero
fuente

para mí coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0")) , así como para coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0",cluster="group"))producir resultados idénticos. ¿Sabes por qué este es el caso?

Peter Pan

1

El enlace people.su.se/~ma/clustering.pdf ya no funciona. ¿Recuerdas el título de la página?

MERose

8

Después de mucha lectura, encontré la solución para hacer clustering dentro del lmmarco.

Hay un excelente libro blanco de Mahmood Arai que proporciona un tutorial sobre la agrupación en clúster en el lmmarco, lo que hace con correcciones de grados de libertad en lugar de mis intentos desordenados anteriores. Él proporciona sus funciones para matrices de covarianza de agrupación de una y dos vías aquí .

Finalmente, aunque el contenido no está disponible de forma gratuita, la Econometría en su mayoría inofensiva de Angrist y Pischke tiene una sección sobre agrupación que fue muy útil.

Actualice el 27/04/2015 para agregar código de la publicación del blog .

api=read.csv("api.csv") #create the variable api from the corresponding csv
attach(api) # attach of data.frame objects
api1=api[c(1:6,8:310),] # one missing entry in row nr. 7
modell.api=lm(API00 ~ GROWTH + EMER + YR_RND, data=api1) # creation of a simple linear model for API00 using the regressors Growth, Emer and Yr_rnd.

##creation of the function according to Arai:
clx <- function(fm, dfcw, cluster) {
    library(sandwich)
    library(lmtest)
    library(zoo)
    M <- length(unique(cluster))
    N <- length(cluster)
    dfc <- (M/(M-1))*((N-1)/(N-fm$rank)) # anpassung der freiheitsgrade
    u <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum))
    vcovCL <-dfc * sandwich (fm, meat = crossprod(u)/N) * dfcw
    coeftest(fm, vcovCL)
}

clx(modell.api, 1, api1$DNUM) #creation of results.

Richard Herron
fuente

1

El periódico de Arai ya no está en línea. ¿Puedes proporcionar el enlace real?

MERose

@MERose - ¡Lo siento! ¡Desafortunadamente no podemos adjuntar archivos PDF! Encontré esta publicación de blog que compara el código. Editaré esta respuesta para incluir el código.

Richard Herron

Esto podría ser una versión actualizada del papel blanco: Mahmood Arai, Cluster-errores estándar robustos utilizando R .

gung - Restablece a Monica

4

La forma más fácil 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 una excelente publicación sobre clustering dentro del marco lm. El sitio también proporciona la función de resumen modificada para la agrupación en clúster de una y dos vías. Puedes encontrar la función y el tutorial aquí .

Tolga Suer
fuente

1

Si no tiene un timeíndice, no necesita uno: plmagregará uno ficticio por sí mismo y no se utilizará a menos que lo solicite. Entonces 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 encontrado un error plm. (Este error ahora se ha solucionado en SVN. Puede instalar la versión de desarrollo desde aquí ).

Pero dado que este sería un timeíndice ficticio 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 plmforma predeterminada, Arellano estimará agrupado por grupos SE. Que es diferente de lo que vcovHC.lm()en sandwichestimará (por ejemplo, los vcovHCPEs en la pregunta original), que es la SE heteroscedasticidad coherente con ninguna agrupación.

Un enfoque separado para esto es apegarse a lmlas regresiones variables ficticias y al paquete multiwayvcov .

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, obtendrá los SE de Arellano (1987) con agrupamiento por grupo. El multiwayvcovpaquete es una evolución directa y significativa de las funciones de agrupamiento originales de Arai.

También puede ver la matriz de varianza-covarianza resultante de ambos enfoques, obteniendo 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

Landroni
fuente

Consulte este enlace: multiwayvcov se deprecia: sites.google.com/site/npgraham1/research/code

HoneyBuddha

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

Respuestas: