Cluster de erro padrão em R (manualmente ou em plm)

33

Estou tentando entender o erro padrão "clustering" e como executar no R (é trivial no Stata). No RI, não obtive sucesso usando plmou escrevendo minha própria função. Vou usar os diamondsdados do ggplot2pacote.

Eu posso fazer efeitos fixos com variáveis fictícias

> 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

ou des-significando os lados esquerdo e direito (sem regressores invariantes no tempo aqui) e corrigindo graus de liberdade.

> # 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

Não consigo replicar esses resultados com plm, porque não tenho um índice de "tempo" (ou seja, esse não é realmente um painel, apenas clusters que poderiam ter um viés comum em termos de erro).

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

Também tentei codificar minha própria matriz de covariância com erro padrão em cluster usando a explicação da Stata sobre sua clusteropção ( explicada aqui ), que é resolver onde , si o número de clusters, é o resíduo para a observação e é o vetor de linha dos preditores, incluindo a constante (isso também aparece como equação (7.22) na seção transversal de Wooldridge e nos dados do painel

{\hat{V}}_{c eu você s t e r} = (X^{'} X)^{- 1} (\sum_{j = 1}^{n_{c}} {você}_{j}^{'} {você}_{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$ ) Mas o código a seguir fornece matrizes de covariância muito grandes. Esses valores são muito altos, devido ao pequeno número de clusters que tenho? Como não consigo plmcriar clusters com um fator, não sei como fazer o benchmark do meu 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

Isso pode ser feito em R? É uma técnica bastante comum em econometria (há um breve tutorial nesta palestra ), mas não consigo descobrir isso em R. Obrigado!

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

7

@ricardh, xingue todos os economistas por não verificarem se o termo que eles querem usar já é usado nas estatísticas. Cluster neste contexto significa grupo e não tem nenhuma relação com a análise de cluster, é por isso que rseek forneceu resultados não relacionados. Por isso, removi a marca de cluster. Para análise dos dados do painel, consulte o pacote plm . Ele tem uma bela vinheta, para que você possa encontrar o que deseja. Quanto à sua pergunta, não está claro o que você deseja. Dentro dos erros padrão do grupo?

precisa saber é o seguinte

@ricardh, ajudaria muito se você pudesse criar um link para algum manual da Stata, onde essa clusteropção é explicada. Estou certo de que seria possível replicar em R.

mpiktas

2

+1 para esse comentário. economistas colonizam a terminologia como loucos. Embora às vezes seja difícil escolher o vilão. Eu demorei um pouco, por exemplo, até perceber factorque não tinha nada a ver com isso, factanalmas se refere a variáveis categorizadas. No entanto, cluster em R se refere à análise de cluster, k-means é apenas o método de particionamento: statmethods.net/advstats/cluster.html . Não entendi sua pergunta, mas também acho que o cluster não tem nada a ver com isso.

hans0l0

@mpiktas, @ ran2 - Obrigado! Espero ter esclarecido a pergunta. Em resumo, por que o "cluster de erro padrão" existe se são apenas efeitos fixos, que já existiam?

Richard Herron

1

A função cluster.vcov no pacote "multiwayvcov" faz o que você está procurando.

29

Para erros padrão brancos agrupados por grupo com a plmestrutura, tente

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

onde model.plmé um modelo plm.

Veja também este link

http://www.inside-r.org/packages/cran/plm/docs/vcovHC ou a documentação do pacote plm

EDITAR:

Para clustering bidirecional (por exemplo, grupo e hora), consulte o seguinte link:

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

Aqui está outro guia útil para o pacote plm especificamente que explica opções diferentes para erros padrão em cluster:

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

Clustering e outras informações, especialmente para Stata, podem ser encontradas aqui:

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

EDIT 2:

Aqui estão exemplos que comparam R e stata: http://www.richard-bluhm.com/clustered-ses-in-r-and-stata-2/

Além disso, multiwayvcovpode ser útil. Esta publicação fornece uma visão geral útil: http://rforpublichealth.blogspot.dk/2014/10/easy-clustered-standard-errors-in-r.html

A partir da documentação:

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)

chandler
fonte

para mim coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0")) , além de coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0",cluster="group"))produzir resultados idênticos. Você sabe por que esse é o caso?

Peter Pan

1

O link people.su.se/~ma/clustering.pdf não está mais funcionando. Você se lembra do título da página?

MERose

8

Depois de muita leitura, encontrei a solução para fazer clustering dentro da lmestrutura.

Há um excelente white paper de Mahmood Arai que fornece um tutorial sobre clustering na lmestrutura, que ele faz com correções de graus de liberdade, em vez de minhas tentativas confusas acima. Ele fornece suas funções para matrizes de covariância de cluster unidirecional e bidirecional aqui .

Finalmente, embora o conteúdo não esteja disponível gratuitamente, o Mostly Harmless Econometrics de Angrist e Pischke possui uma seção sobre clustering que foi muito útil.

Atualize em 27/04/2015 para adicionar código da postagem do 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
fonte

1

O artigo de Arai não está mais online. Você pode fornecer o link real?

MERose

@MERose - Desculpe! Infelizmente não podemos anexar PDFs! Encontrei este post de blog que compara o código. Vou editar esta resposta para incluir o código.

Richard Herron

Esta poderia ser uma versão atualizada do Livro Branco: Mahmood Arai, erros padrão Cluster-robustos usando R .

gung - Restabelece Monica

4

A maneira mais fácil de calcular erros padrão em cluster no R é usar a função de resumo modificada.

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

Há um excelente post sobre clustering na estrutura lm. O site também fornece a função de resumo modificada para clustering unidirecional e bidirecional. Você pode encontrar a função e o tutorial aqui .

Tolga Suer
fonte

1

Se você não possui um timeíndice, não precisa de um: plmadicionará um fictício por si só e não será usado a menos que você o solicite. Portanto, esta chamada deve funcionar :

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

Exceto que não, o que sugere que você encontrou um bug plm. (Este bug foi corrigido no SVN. Você pode instalar a versão de desenvolvimento a partir daqui .)

Mas como esse seria um timeíndice fictício de qualquer maneira, podemos criá-lo por nós mesmos:

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

Agora isso 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()na plmvontade pela estimativa padrão Arellano agrupados por SEs grupo. Que é diferente do que vcovHC.lm()em sandwichestimará (por exemplo, as vcovHCSEs na pergunta original), que é SEs heterocedasticidade consistente sem clustering.

Uma abordagem separada disso é a lmregressão de variáveis fictícias e o pacote 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

Nos dois casos, você obterá os SEs de Arellano (1987) com agrupamento por grupo. O multiwayvcovpacote é uma evolução direta e significativa das funções de clustering originais do Arai.

Você também pode examinar a matriz de variância-covariância resultante de ambas as abordagens, produzindo a mesma estimativa de variância 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
fonte

Por favor, veja este link: multiwayvcov está depreciado: sites.google.com/site/npgraham1/research/code

HoneyBuddha

Cluster de erro padrão em R (manualmente ou em plm)

Respostas: