Standardowe grupowanie błędów w R (ręcznie lub w trybie plm)

33

Próbuję zrozumieć standardowy błąd „klastrowanie” i sposób wykonania w języku R (w Stacie jest to trywialne). W RI nie udało mi się ani użyć ani plmnapisać własnej funkcji. Użyję diamondsdanych z ggplot2paczki.

Potrafię robić stałe efekty z dowolnymi zmiennymi obojętnymi

> 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 

lub przez pozbawienie znaczenia zarówno lewej, jak i prawej strony (brak regresorów niezmiennych w czasie) i korygowanie stopni swobody.

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

Nie mogę powielić tych wyników plm, ponieważ nie mam indeksu „czasu” (tzn. To nie jest tak naprawdę panel, tylko klastry, które mogą mieć wspólne odchylenie pod względem błędów).

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

Próbowałem także zakodować własną macierz kowariancji za pomocą standardowego błędu klastrowego, używając wyjaśnień Staty dotyczących ich clusteropcji ( wyjaśnionej tutaj ), która polega na rozwiązaniu gdzie , si liczba klastrów, jest resztą dla obserwacji a jest wektorem rzędów predyktorów, w tym stałą (pojawia się to również jako równanie (7.22) w przekroju Wooldridge'a i danych panelowych

V^cluster=(XX)1(j=1ncujuj)(XX)1
uj=cluster jeixinceiithxi). Ale poniższy kod podaje bardzo duże macierze kowariancji. Czy te bardzo duże wartości są biorąc pod uwagę małą liczbę klastrów, które mam? Biorąc pod uwagę, że nie mogę tworzyć plmklastrów na jednym czynniku, nie jestem pewien, jak przeprowadzić test porównawczy mojego kodu.
> # 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

Czy można to zrobić w R? Jest to dość powszechna technika w ekonometrii (w tym wykładzie znajduje się krótki samouczek ), ale nie mogę jej zrozumieć w R. Dzięki!

Richard Herron
źródło
7
@ricardh, przeklinaj wszystkich ekonomistów, że nie sprawdzili, czy termin, którego chcą użyć, jest już używany w statystykach. Klaster w tym kontekście oznacza grupę i jest całkowicie niezwiązany z analizą skupień, dlatego rseek dał ci niepowiązane wyniki. Dlatego usunąłem tag grupowania. Do analizy danych panelowych sprawdź pakiet plm . Ma ładną winietę, więc możesz znaleźć to, czego chcesz. Jeśli chodzi o twoje pytanie, nie jest jasne, czego chcesz. W obrębie standardowych błędów grupy?
mpiktas
@ricardh, byłoby bardzo pomocne, gdybyś mógł link do instrukcji Staty, w której ta clusteropcja jest wyjaśniona. Jestem pewien, że replikacja w R. byłaby możliwa
mpiktas,
2
+1 za ten komentarz. ekonomiści kolonizują terminologię jak szaleni. Chociaż czasami trudno jest wybrać złoczyńcę. Poświęciłem trochę czasu, np. Dopóki nie zdałem sobie sprawy, factorże nie mam nic wspólnego, factanalale odnosi się do zmiennych skategoryzowanych. Jakkolwiek klaster w R odnosi się do analizy skupień, k-średnich jest tylko metodą partycjonowania: statmethods.net/advstats/cluster.html . Nie dostaję twojego pytania, ale wydaje mi się, że klaster nie ma z tym nic wspólnego.
hans0l0
@mpiktas, @ ran2 - Dzięki! Mam nadzieję, że wyjaśniłem pytanie. Krótko mówiąc, dlaczego istnieje „standardowe grupowanie błędów”, jeśli są to tylko ustalone efekty, które już istniały?
Richard Herron
1
Funkcjauster.vcov w pakiecie „multiwayvcov” wykonuje to, czego szukasz.

Odpowiedzi:

29

W przypadku białych standardowych błędów skupionych w grupach według plmstruktury spróbuj

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

gdzie model.plmjest model plm.

Zobacz także ten link

http://www.inside-r.org/packages/cran/plm/docs/vcovHC lub dokumentacja pakietu plm

EDYTOWAĆ:

W przypadku klastrowania dwukierunkowego (np. Grupa i czas) zobacz następujący link:

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

Oto kolejny pomocny przewodnik dla pakietu PLM, który wyjaśnia różne opcje standardowych błędów klastrowanych:

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

Klastry i inne informacje, szczególnie dla Staty, można znaleźć tutaj:

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

EDYCJA 2:

Oto przykłady, które porównują R i stata: http://www.richard-bluhm.com/clustered-ses-in-r-and-stata-2/

Ponadto multiwayvcovmogą być pomocne. Ten post zawiera pomocne omówienie: http://rforpublichealth.blogspot.dk/2014/10/easy-clustered-standard-errors-in-r.html

Z dokumentacji:

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)
kupiec
źródło
dla mnie, coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0")) a także coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0",cluster="group"))produkować identyczne wyniki. Czy wiesz, dlaczego tak jest?
Peter Pan,
1
Link people.su.se/~ma/clustering.pdf nie działa. Czy pamiętasz tytuł strony?
MERose
8

Po wielu lekturach znalazłem rozwiązanie do tworzenia klastrów w lmramach.

Istnieje doskonała biała księga Mahmooda Arai, która zawiera samouczek na temat grupowania w lmramach, który wykonuje z korektami stopni swobody zamiast moich niechlujnych prób powyżej. On zapewnia swoje funkcje zarówno jedno- i dwukierunkowych macierzy kowariancji klastrów tutaj .

Wreszcie, chociaż treść nie jest dostępna za darmo, Najbardziej nieszkodliwe ekonometria” autorstwa Angrista i Pischke ma bardzo pomocny rozdział dotyczący grupowania.


Zaktualizuj 27.04.2015, aby dodać kod z posta na blogu .

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
źródło
1
Artykuł Arai nie jest już dostępny w Internecie. Czy możesz podać rzeczywisty link?
MERose
@MERose - Przepraszamy! Niestety nie możemy załączyć plików PDF! Znalazłem ten post na blogu, który porównuje kod. Zmienię tę odpowiedź, aby dołączyć kod.
Richard Herron
To może być zaktualizowana wersja białego papieru: Mahmood Arai, Cluster-solidne błędy standardowe wykorzystujące R .
gung - Przywróć Monikę
4

Najłatwiejszym sposobem obliczenia klastrowanych błędów standardowych w R jest użycie zmodyfikowanej funkcji podsumowania.

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

Istnieje doskonały post na temat klastrowania w ramach lm. Witryna udostępnia również zmodyfikowaną funkcję podsumowania dla grupowania jedno- i dwukierunkowego. Możesz znaleźć funkcję i samouczek tutaj .

Tolga Suer
źródło
1

Jeśli nie masz timeindeksu, nie potrzebujesz go: sam plmdoda fikcyjny i nie będzie używany, dopóki go nie poprosisz. To połączenie powinno działać :

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

Tyle że nie, co sugeruje, że trafiłeś w błąd plm. (Ten błąd został już naprawiony w SVN. Możesz zainstalować wersję rozwojową tutaj .)

Ponieważ jednak i tak byłby to fikcyjny timeindeks, możemy go stworzyć sami:

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

Teraz to działa:

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

Ważna uwaga : vcovHC.plm()w plmwoli domyślnie szacunków Arellano skupione przez SE grupowych. Który jest inny od tego, co vcovHC.lm()w sandwichoszacuje (np vcovHCSE w oryginalnym pytanie), który jest zgodny Heteroskedastyczność-SE bez grupowania.


Odrębnym podejściem do tego jest trzymanie się lmfałszywych regresji zmiennych i pakietu 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

W obu przypadkach otrzymasz SE Arellano (1987) z grupowaniem według grup. multiwayvcovPakiet jest bezpośredni i znaczący rozwój pierwotnych funkcji grupowania Arai jest.

Możesz także spojrzeć na wynikową macierz wariancji-kowariancji z obu podejść, uzyskując takie same oszacowanie wariancji dla 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
źródło
Proszę zobaczyć ten link: multiwayvcov jest amortyzowany: sites.google.com/site/npgraham1/research/code
HoneyBuddha