Стандартная кластеризация ошибок в R (вручную или в plm)

33

Я пытаюсь понять стандартную ошибку «кластеризация» и как выполнить в R (это тривиально в Stata). В РИ были неудачные попытки использования либо plmнаписания моей собственной функции. Я буду использовать diamondsданные из ggplot2пакета.

Я могу сделать фиксированные эффекты с помощью фиктивных переменных

> 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

или путем де-значения как левой, так и правой частей (здесь нет времени-инвариантных регрессоров) и корректировки степеней свободы.

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

Я не могу воспроизвести эти результаты plm, потому что у меня нет индекса «времени» (т. Е. На самом деле это не панель, а кластеры, которые могут иметь общее смещение в терминах ошибок).

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

Я также попытался закодировать свою собственную ковариационную матрицу с кластеризованной стандартной ошибкой, используя объяснение Stata их clusterопции ( объясняется здесь ), которая должна решить где , si количество кластеров, - остаток для наблюдения а - вектор строки предикторов, включая константу (это также появляется в виде уравнения (7.22) в поперечном сечении Вулдриджа и панельных данных

{\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$ ). Но следующий код дает очень большие ковариационные матрицы. Являются ли эти очень большие значения с учетом небольшого количества кластеров, которые у меня есть? Учитывая, что я не могу plmсделать кластеры по одному фактору, я не уверен, как сравнить мой код.

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

Это можно сделать в R? Это довольно распространенный метод в эконометрике (в этой лекции есть краткое руководство ), но я не могу понять это в R. Спасибо!

r panel-data standard-error fixed-effects-model clustered-standard-errors Ричард Херрон
источник

7

@ricardh, прокляни всех экономистов за то, что они не проверяют, используется ли уже используемый ими термин в статистике. Кластер в этом контексте означает группу и совершенно не связан с кластерным анализом, поэтому rseek дал вам не связанные результаты. Поэтому я удалил тег кластеризации. Для анализа данных панели проверьте пакет plm . У него хорошая виньетка, так что вы можете найти то, что хотите. Что касается твоего вопроса, то не понятно, чего ты хочешь. Внутри группы стандартные ошибки?

mpiktas

@ricardh, было бы очень полезно, если бы вы могли сослаться на какое-то руководство по Stata, где clusterобъясняется эта опция. Я уверен, что можно было бы повторить в R.

mpiktas

2

+1 за этот комментарий. экономисты колонизируют терминологию как сумасшедшую. Хотя иногда трудно выбрать злодея. Например, пока я не понял, factorчто это не имеет ничего общего, factanalа относится к категорированным переменным. Однако кластер в R относится к кластерному анализу, а k-means - это просто метод разделения: statmethods.net/advstats/cluster.html . Я не понимаю ваш вопрос, но я также предполагаю, что кластер не имеет к этому никакого отношения.

hans0l0

@mpiktas, @ ran2 - Спасибо! Надеюсь, я прояснил вопрос. Короче говоря, почему существует «стандартная кластеризация ошибок», если это просто фиксированные эффекты, которые уже существовали?

Ричард Херрон

1

Функция cluster.vcov в пакете «multiwayvcov» делает то, что вы ищете.

29

Для белых стандартных ошибок, сгруппированных по группам со plmструктурой, попробуйте

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

где model.plmмодель PLM.

Смотрите также эту ссылку

http://www.inside-r.org/packages/cran/plm/docs/vcovHC или документация пакета plm

РЕДАКТИРОВАТЬ:

Для двусторонней кластеризации (например, группа и время) см. Следующую ссылку:

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

Вот еще одно полезное руководство для пакета plm, в котором объясняются различные варианты кластерных стандартных ошибок:

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

Информация о кластерах и другая информация, особенно для Stata, можно найти здесь:

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

РЕДАКТИРОВАТЬ 2:

Вот примеры, которые сравнивают R и stata: http://www.richard-bluhm.com/clustered-ses-in-r-and-stata-2/

Кроме того, multiwayvcovможет быть полезным. Этот пост содержит полезный обзор: http://rforpublichealth.blogspot.dk/2014/10/easy-clustered-standard-errors-in-r.html

Из документации:

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)

лавочник
источник

для меня, coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0")) а также coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0",cluster="group"))дают идентичные результаты. Вы знаете, почему это так?

Питер Пэн

1

Ссылка people.su.se/~ma/clustering.pdf больше не работает. Вы помните заголовок страницы?

MERose

8

После долгих чтений я нашел решение для кластеризации в lmрамках.

Есть превосходный технический документ Махмуда Араи, который предоставляет учебник по кластеризации в lmрамках, который он делает с поправками на степени свободы вместо моих грязных попыток выше. Он предоставляет свои функции как для одно- и двухсторонних матриц ковариации кластеризации здесь .

Наконец, хотя контент не доступен бесплатно, Angrist и Pischke в Mostly Harmless Econometrics есть раздел о кластеризации, который был очень полезен.

Обновление от 27.04.2015 для добавления кода из поста в блоге .

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.

Ричард Херрон
источник

1

Бумаги Араи больше нет в сети. Можете ли вы предоставить актуальную ссылку?

MERose

@Mose - прости! К сожалению, мы не можем прикрепить PDF-файлы! Я нашел этот пост в блоге, который тестирует код. Я отредактирую этот ответ, чтобы включить код.

Ричард Херрон

Это может быть обновленная версия белой бумаги: Махмуд Arai, Кластер надежные стандартные ошибки с использованием R .

gung - Восстановить Монику

4

Самый простой способ вычислить кластерные стандартные ошибки в R - это использовать модифицированную функцию суммирования.

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

Есть отличный пост по кластеризации в рамках lm. На сайте также представлена измененная функция сводки как для односторонней, так и для двусторонней кластеризации. Вы можете найти функцию и учебник здесь .

Толга Суер
источник

1

Если у вас нет timeиндекса, он вам не нужен: plmон сам добавит фиктивный, и он не будет использоваться, если вы не попросите его. Так что этот вызов должен работать :

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

За исключением того, что это не так, что предполагает, что вы попали в ошибку plm. (Эта ошибка теперь исправлена в SVN. Вы можете установить версию для разработки отсюда .)

Но поскольку в timeлюбом случае это фиктивный индекс, мы можем создать его сами:

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

Теперь это работает:

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

Важное примечание : vcovHC.plm()в plmпо умолчанию будет оценкой Арельяно кластерного групповой СЭП. Который отличается от того, что vcovHC.lm()в sandwichоценит (например, vcovHCпоисковикам в исходном вопросе), то есть гетероскедастичности-последовательной СЭП с не кластеризация.

Отдельный подход к этому - придерживаться lmфиктивных переменных регрессий и пакета 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

В обоих случаях вы получите SE Arellano (1987) с кластеризацией по группам. multiwayvcovПакет является прямой и значительной эволюцией оригинальных функций кластеризации Араи.

Вы также можете посмотреть на полученную матрицу дисперсии-ковариации из обоих подходов, получая одинаковую оценку дисперсии для 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
источник

Пожалуйста, посмотрите эту ссылку: multiwayvcov

устарела

Стандартная кластеризация ошибок в R (вручную или в plm)

Ответы: