Standardfehler-Clustering in R (entweder manuell oder in plm)

33

Ich versuche zu verstehen, Standardfehler "Clustering" und wie in R auszuführen (es ist in Stata trivial). In RI ist es mir nicht gelungen, entweder plmmeine eigene Funktion zu verwenden oder zu schreiben. Ich werde die diamondsDaten aus dem ggplot2Paket verwenden.

Ich kann feste Effekte mit beiden Dummy-Variablen machen

> 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

oder indem Sie die Bedeutung sowohl der linken als auch der rechten Seite aufheben (keine zeitinvarianten Regressoren hier) und Freiheitsgrade korrigieren.

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

Ich kann diese Ergebnisse nicht mit replizieren plm, da ich keinen "Zeit" -Index habe (dh dies ist nicht wirklich ein Panel, sondern nur Cluster, die in ihren Fehlerbegriffen eine gemeinsame Tendenz aufweisen könnten).

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

Ich habe auch versucht, meine eigene Kovarianzmatrix mit geclusterten Standardfehlern zu codieren, indem ich die Erklärung von Stata zu ihrer clusterOption ( hier erklärt ) verwendet habe: wobei , si die Anzahl von Clustern ist, der Rest ist für die Beobachtung und ist der Zeilenvektor von Prädiktoren, einschließlich der Konstanten (dies erscheint auch als Gleichung (7.22) in Wooldridges Querschnitts- und Paneldaten

{\hat{V}}_{c l u s t e r} = (X^{'} X)^{- 1} (\sum_{j = 1}^{n_{c}} u_{j}^{'} u_{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$ ). Der folgende Code liefert jedoch sehr große Kovarianzmatrizen. Sind diese sehr großen Werte angesichts der geringen Anzahl von Clustern, die ich habe? Da ich nicht plmin der Lage bin, Cluster zu einem bestimmten Faktor zu erstellen, bin ich mir nicht sicher, wie ich meinen Code bewerten soll.

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

Kann das in R gemacht werden? Es ist eine in der Ökonometrie weit verbreitete Technik (in dieser Vorlesung gibt es ein kurzes Tutorial ), aber ich kann es in R nicht herausfinden. Danke!

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

7

@ricardh, verfluche alle Ökonomen, die nicht überprüft haben, ob der Begriff, den sie verwenden möchten, bereits in der Statistik verwendet wird. Cluster bedeutet in diesem Zusammenhang Gruppe und steht in keiner Beziehung zur Clusteranalyse. Aus diesem Grund hat rseek Ihnen Ergebnisse geliefert, die in keiner Beziehung stehen. Daher habe ich das Clustering-Tag entfernt. Für die Analyse der Paneldaten siehe Paket plm . Es hat eine schöne Vignette, so dass Sie vielleicht finden, was Sie wollen. Bei Ihrer Frage ist nicht klar, was Sie wollen. Innerhalb der Gruppe Standardfehler?

mpiktas

@ricardh, es wäre sehr hilfreich, wenn Sie auf ein Handbuch von Stata verweisen könnten, in dem diese clusterOption erklärt wird. Ich bin sicher, es wäre möglich, in R.

mpiktas

2

+1 für diesen Kommentar. Ökonomen kolonisieren die Terminologie wie verrückt. Obwohl es manchmal schwierig ist, den Bösewicht zu finden. Ich factorhabe eine Weile gebraucht, bis mir klar wurde, dass ich nichts damit zu tun hatte, factanalsondern mich auf kategorisierte Variablen bezog. Cluster in R bezieht sich jedoch auf die Clusteranalyse, k-means ist nur DIE Partitionierungsmethode: statmethods.net/advstats/cluster.html . Ich verstehe Ihre Frage nicht, aber ich schätze, Cluster hat nichts damit zu tun.

hans0l0

@mpiktas, @ ran2 - Danke! Ich hoffe, ich habe die Frage geklärt. Kurz gesagt, warum gibt es "Standard Error Clustering", wenn es sich nur um festgelegte Effekte handelt, die es bereits gab?

Richard Herron

1

Die cluster.vcov-Funktion im "multiwayvcov" -Paket macht das, wonach Sie suchen.

29

Bei White-Standardfehlern gruppiert nach Gruppe mit dem plmFramework versuchen

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

Wo model.plmist ein plm-Modell.

Siehe auch diesen Link

http://www.inside-r.org/packages/cran/plm/docs/vcovHC oder die plm-Paketdokumentation

BEARBEITEN:

Informationen zum bidirektionalen Clustering (z. B. Gruppe und Zeit) finden Sie unter folgendem Link:

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

Hier ist eine weitere hilfreiche Anleitung für das plm-Paket, in der verschiedene Optionen für Clustered-Standardfehler erläutert werden:

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

Clustering und andere Informationen, insbesondere zu Stata, finden Sie hier:

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

EDIT 2:

Hier sind Beispiele, die R und Stata vergleichen: http://www.richard-bluhm.com/clustered-ses-in-r-and-stata-2/

Auch das multiwayvcovkann hilfreich sein. Dieser Beitrag bietet eine hilfreiche Übersicht: http://rforpublichealth.blogspot.dk/2014/10/easy-clustered-standard-errors-in-r.html

Aus der Dokumentation:

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)

Ausrüster
quelle

für mich coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0")) genauso wie coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0",cluster="group"))für identische ergebnisse. Wissen Sie, warum das so ist?

Peter Pan

1

Der Link people.su.se/~ma/clustering.pdf funktioniert nicht mehr. Erinnerst du dich an den Titel der Seite?

MERose

8

Nach vielem Lesen fand ich die Lösung für das Clustering innerhalb des lmFrameworks.

Es gibt ein exzellentes Whitepaper von Mahmood Arai, das ein Tutorial zum Clustering im lmFramework enthält, das er mit Korrekturen der Freiheitsgrade anstelle meiner oben beschriebenen chaotischen Versuche durchführt. Er stellt seine Funktionen sowohl für Ein- und Zwei-Wege - Clustering Kovarianzmatrizen hier .

Obwohl der Inhalt nicht kostenlos verfügbar ist, hat Angrist und Pischkes Mostly Harmless Econometrics einen Abschnitt zum Clustering, der sehr hilfreich war.

Update am 27.04.2015, um Code aus dem Blogeintrag hinzuzufügen .

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
quelle

1

Arais Artikel ist nicht mehr online. Können Sie den tatsächlichen Link angeben?

MERose

@MERose - Entschuldigung! Leider können wir keine PDFs anhängen! Ich habe diesen Blog-Beitrag gefunden , der den Code bewertet. Ich werde diese Antwort so bearbeiten, dass sie den Code enthält.

Richard Herron

: Dies könnte eine aktualisierte Version des Weißbuch sein Mahmood Arai, Cluster-robuster Standardfehler unter Verwendung von R .

gung - Wiedereinsetzung von Monica

4

Der einfachste Weg, gruppierte Standardfehler in R zu berechnen, ist die Verwendung der geänderten Zusammenfassungsfunktion.

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

Es gibt einen ausgezeichneten Beitrag zum Clustering im Rahmen von lm. Die Site bietet auch die geänderte Zusammenfassungsfunktion für Einweg- und Zweiweg-Clustering. Die Funktion und das Tutorial finden Sie hier .

Tolga Suer
quelle

1

Wenn Sie keinen timeIndex haben, brauchen Sie keinen: Er plmfügt selbst einen fiktiven hinzu und wird nur verwendet, wenn Sie danach fragen. Also sollte dieser Aufruf funktionieren :

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

Nur dass dies nicht der Fall ist, was darauf hindeutet, dass Sie einen Fehler gefunden haben plm. (Dieser Fehler wurde nun in SVN behoben. Sie können die Entwicklungsversion von hier aus installieren .)

Da dies jedoch ohnehin ein fiktiver timeIndex wäre, können wir ihn selbst erstellen:

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

Das funktioniert jetzt:

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

Wichtiger Hinweis : vcovHC.plm()in plmwird standardmäßig Schätzung Arellano gruppierten von Gruppe SEs. Das ist verschieden von dem, was vcovHC.lm()in sandwichschätzen werden (zB die vcovHCSEs in der ursprünglichen Frage), das heteroscedasticity konsistenten SEs ohne Clustering ist.

Ein separater Ansatz besteht darin, an lmDummy-Variablen-Regressionen und dem multiwayvcov- Paket festzuhalten .

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

In beiden Fällen erhalten Sie die Arellano (1987) SEs mit Clustering nach Gruppe. Das multiwayvcovPaket ist eine direkte und signifikante Weiterentwicklung der ursprünglichen Clustering-Funktionen von Arai.

Sie können auch die resultierende Varianz-Kovarianz-Matrix aus beiden Ansätzen betrachten und die gleiche Varianzschätzung erhalten für 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
quelle

Bitte sehen Sie diesen Link: multiwayvcov wird abgeschrieben: sites.google.com/site/npgraham1/research/code

HoneyBuddha

Standardfehler-Clustering in R (entweder manuell oder in plm)

Antworten: