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

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

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

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)

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.

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 .

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