Pengelompokan kesalahan standar dalam R (baik secara manual atau dalam plm)

33

Saya mencoba memahami kesalahan "clustering" standar dan bagaimana mengeksekusi di R (itu sepele di Stata). Di RI telah gagal menggunakan salah satu plmatau menulis fungsi saya sendiri. Saya akan menggunakan diamondsdata dari ggplot2paket.

Saya bisa melakukan efek tetap dengan salah satu variabel dummy

> 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

atau dengan mengartikan kedua sisi kiri dan kanan (tidak ada waktu yang regresi di sini) dan memperbaiki derajat kebebasan.

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

Saya tidak dapat mereplikasi hasil ini dengan plm, karena saya tidak memiliki indeks "waktu" (yaitu, ini bukan benar-benar panel, hanya cluster yang bisa memiliki bias umum dalam hal kesalahan mereka).

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

Saya juga mencoba untuk kode matriks kovarians saya sendiri dengan kesalahan standar berkerumun menggunakan penjelasan Stata tentang clusteropsi mereka ( dijelaskan di sini ), yaitu untuk memecahkan mana , si jumlah cluster, adalah residual untuk observasi dan adalah vektor baris prediktor, termasuk konstanta (ini juga muncul sebagai persamaan (7.22) dalam Cross Section dan Panel Data Wooldridge

{\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$ ). Tetapi kode berikut memberikan matriks kovarian yang sangat besar. Apakah nilai-nilai yang sangat besar ini diberikan sejumlah kecil cluster yang saya miliki? Mengingat bahwa saya tidak plmdapat melakukan cluster pada satu faktor, saya tidak yakin bagaimana melakukan benchmark kode saya.

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

Bisakah ini dilakukan di R? Ini adalah teknik yang cukup umum dalam ekonometrik (ada tutorial singkat dalam kuliah ini ), tapi saya tidak bisa mengetahuinya dalam R. Terima kasih!

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

7

@ricardh, kutuk semua ekonom karena tidak memeriksa apakah istilah yang ingin mereka gunakan sudah digunakan dalam statistik. Cluster dalam konteks ini berarti grup dan sama sekali tidak terkait dengan analisis cluster, inilah mengapa rseek memberi Anda hasil yang tidak terkait. Karenanya saya menghapus tag pengelompokan. Untuk analisis data panel, periksa plm paket . Ini memiliki sketsa yang bagus, sehingga Anda dapat menemukan apa yang Anda inginkan. Adapun pertanyaan Anda tidak jelas apa yang Anda inginkan. Dalam kesalahan standar grup?

mpiktas

@ricardh, akan sangat membantu jika Anda dapat menautkan ke beberapa manual Stata di mana clusteropsi ini dijelaskan. Saya yakin akan mungkin untuk mereplikasi di R.

mpiktas

2

+1 untuk komentar itu. ekonom mengolonisasi terminologi seperti orang gila. Meskipun terkadang sulit untuk memilih penjahatnya. Ii mengambil beberapa saat misalnya sampai saya menyadari factortidak ada hubungannya dengan factanaltetapi mengacu pada variabel yang dikategorikan. Namun cluster dalam R mengacu pada analisis cluster, k-means hanyalah metode partisi: statmethods.net/advstats/cluster.html . Saya tidak mendapatkan pertanyaan Anda, tapi saya kira cluster tidak ada hubungannya dengan itu.

hans0l0

@mpiktas, @ ran2 - Terima kasih! Saya harap saya mengklarifikasi pertanyaan itu. Singkatnya, mengapa "pengelompokan kesalahan standar" ada jika hanya efek tetap, yang sudah ada?

Richard Herron

1

Fungsi cluster.vcov dalam paket "multiwayvcov" melakukan apa yang Anda cari.

29

Untuk kesalahan standar putih dikelompokkan berdasarkan grup dengan plmkerangka coba

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

di mana model.plmmodel PLM.

Lihat juga tautan ini

http://www.inside-r.org/packages/cran/plm/docs/vcovHC atau dokumentasi paket plm

EDIT:

Untuk pengelompokan dua arah (mis. Grup dan waktu) lihat tautan berikut:

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

Berikut ini adalah panduan bermanfaat lain untuk paket PLM yang secara khusus menjelaskan berbagai opsi untuk kesalahan standar berkerumun:

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

Pengelompokan dan informasi lainnya, terutama untuk Stata, dapat ditemukan di sini:

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

EDIT 2:

Berikut adalah contoh yang membandingkan R dan stata: http://www.richard-bluhm.com/clustered-ses-in-r-and-stata-2/

Juga, multiwayvcovmungkin bermanfaat. Posting ini memberikan ikhtisar bermanfaat: http://rforpublichealth.blogspot.dk/2014/10/easy-clustered-standard-errors-in-r.html

Dari dokumentasi:

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)

pedagang lilin
sumber

bagi saya coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0")) juga coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0",cluster="group"))menghasilkan hasil yang identik. Apakah Anda tahu mengapa demikian?

Peter Pan

1

Tautan people.su.se/~ma/clustering.pdf tidak lagi berfungsi. Apakah Anda ingat judul halaman?

MERose

8

Setelah banyak membaca, saya menemukan solusi untuk melakukan pengelompokan dalam lmkerangka kerja.

Ada buku putih yang sangat bagus dari Mahmood Arai yang menyediakan tutorial tentang pengelompokan dalam lmkerangka kerja, yang dia lakukan dengan koreksi derajat kebebasan alih-alih upaya saya yang berantakan di atas. Ia menyediakan fungsinya untuk matriks kovarians pengelompokan satu dan dua arah di sini .

Akhirnya, meskipun konten tidak tersedia gratis, Angrist dan Mostly Harmless Econometrics dari Pischke memiliki bagian tentang pengelompokan yang sangat membantu.

Perbarui pada 4/27/2015 untuk menambahkan kode dari posting 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
sumber

1

Koran Arai tidak lagi online. Bisakah Anda memberikan tautan yang sebenarnya?

MERose

@MERose - Maaf! Sayangnya kami tidak dapat melampirkan pdf! Saya menemukan posting blog ini yang menandai kode. Saya akan mengedit jawaban ini untuk memasukkan kode.

Richard Herron

Ini mungkin versi terbaru dari kertas putih: Mahmood Arai, kesalahan standar Cluster-kuat menggunakan R .

gung - Reinstate Monica

4

Cara termudah untuk menghitung kesalahan standar berkerumun di R adalah dengan menggunakan fungsi ringkasan yang dimodifikasi.

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

Ada posting yang sangat baik tentang pengelompokan dalam kerangka lm. Situs ini juga menyediakan fungsi ringkasan yang dimodifikasi untuk pengelompokan satu dan dua arah. Anda dapat menemukan fungsi dan tutorialnya di sini .

Tolga Suer
sumber

1

Jika Anda tidak memiliki timeindeks, Anda tidak perlu indeks: plmdengan sendirinya akan menambahkan indeks fiktif, dan itu tidak akan digunakan kecuali Anda memintanya. Jadi panggilan ini harus berfungsi :

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

Kecuali itu tidak, yang berarti Anda telah menemukan bug plm. (Bug ini sekarang telah diperbaiki di SVN. Anda dapat menginstal versi pengembangan dari sini .)

Tetapi karena ini akan menjadi timeindeks fiktif , kita dapat membuatnya sendiri:

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

Sekarang ini berfungsi:

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

Catatan penting : vcovHC.plm()dalam plmperkiraan default Arellano akan dikelompokkan menurut SE kelompok. Yang berbeda dari apa yang vcovHC.lm()di sandwichestimasi akan (misalnya vcovHCSE dalam pertanyaan awal), yaitu SE heteroscedasticity-konsisten tanpa clustering.

Pendekatan terpisah untuk ini adalah menempel pada lmregresi variabel dummy dan paket 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

Dalam kedua kasus Anda akan mendapatkan Arellano (1987) SE dengan pengelompokan berdasarkan kelompok. The multiwayvcovpaket evolusi langsung dan signifikan dari fungsi pengelompokan asli Arai.

Anda juga dapat melihat matriks varians-kovarians yang dihasilkan dari kedua pendekatan, menghasilkan estimasi varians yang sama untuk 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
sumber

Silakan lihat tautan ini: multiwayvcov disusutkan: sites.google.com/site/npgraham1/research/code

HoneyBuddha

Pengelompokan kesalahan standar dalam R (baik secara manual atau dalam plm)

Jawaban: