Regroupement des erreurs standard dans R (manuellement ou dans plm)

33

J'essaie de comprendre le "clustering" d'erreur standard et comment exécuter dans R (c'est trivial dans Stata). En RI ont été infructueux en utilisant plmou en écrivant ma propre fonction. Je vais utiliser les diamondsdonnées du ggplot2paquet.

Je peux faire des effets fixes avec des variables factices

> 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

ou en dé-signifiant les côtés gauche et droit (pas de régresseurs invariants dans le temps ici) et en corrigeant les degrés de liberté.

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

Je ne peux pas reproduire ces résultats avec plm, car je n'ai pas d'indice "temps" (c'est-à-dire, ce n'est pas vraiment un panel, mais seulement des grappes qui pourraient avoir un biais commun dans leurs termes d'erreur).

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

J'ai aussi essayé de coder ma propre matrice de covariance avec une erreur standard groupée en utilisant l'explication de Stata de leur clusteroption ( expliquée ici ), qui consiste à résoudre où , si le nombre de clusters, est le résidu pour l' observation et est le vecteur de rangée des prédicteurs, y compris la constante (elle apparaît également sous la forme de l'équation (7.22) dans les données de la section transversale et du panel de 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$ ). Mais le code suivant donne de très grandes matrices de covariance. Ces valeurs sont-elles très grandes compte tenu du petit nombre de grappes que j'ai? Étant donné que je ne parviens pas plmà créer des clusters sur un facteur, je ne sais pas comment analyser mon code.

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

Cela peut-il être fait en R? C'est une technique assez courante en économétrie (il y a un bref tutoriel dans cette leçon ), mais je ne peux pas le comprendre en R. Merci!

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

7

@ricardh, insulter tous les économistes pour ne pas avoir vérifié si le terme qu'ils veulent utiliser est déjà utilisé dans les statistiques. Cluster dans ce contexte signifie groupe et n'a aucun lien avec l'analyse en cluster, c'est pourquoi rseek vous a donné des résultats sans lien. Par conséquent, j'ai supprimé la balise de clustering. Pour l'analyse des données de panel, consultez le package plm . Il a une belle vignette, de sorte que vous pouvez trouver ce que vous voulez. Quant à votre question, ce que vous voulez n'est pas clair. Dans les erreurs standard de groupe?

Mpiktas

@ricardh, cela aiderait beaucoup si vous pouviez créer un lien vers un manuel de Stata où cette clusteroption est expliquée. Je suis sûr qu'il serait possible de reproduire R.

mpiktas le

2

+1 pour ce commentaire. les économistes colonisent la terminologie comme un fou. Bien qu'il soit parfois difficile de choisir le méchant. J'ai pris un certain temps, par exemple, jusqu'à ce que je réalise que je factorn'avais rien à voir avec des factanalvariables catégorisées. Cependant, cluster dans R fait référence à l'analyse de cluster, k-means est simplement LA méthode de partitionnement: statmethods.net/advstats/cluster.html . Je ne comprends pas votre question, mais je suppose aussi que le cluster n’a rien à voir avec cela.

hans0l0

@mpiktas, @ ran2 - Merci! J'espère avoir clarifié la question. En bref, pourquoi existe-t-il un "groupement d’erreurs standard" s’il s’agit uniquement d’effets fixes qui existaient déjà?

Richard Herron

1

La fonction cluster.vcov du package "multiwayvcov" fait ce que vous recherchez.

29

Pour les erreurs standard blanches regroupées par groupe avec le plmframework, essayez

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

où model.plmest un modèle plm.

Voir aussi ce lien

http://www.inside-r.org/packages/cran/plm/docs/vcovHC ou la documentation du package plm

MODIFIER:

Pour un regroupement bidirectionnel (par exemple, groupe et heure), voir le lien suivant:

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

Voici un autre guide utile pour le package plm qui explique différentes options pour les erreurs standard en cluster:

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

Le regroupement et d'autres informations, en particulier pour Stata, peuvent être trouvés ici:

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

EDIT 2:

Voici des exemples qui comparent R et stata: http://www.richard-bluhm.com/clustered-ses-in-r-and-stata-2/

En outre, le multiwayvcovpeut être utile. Cet article fournit un aperçu utile: http://rforpublichealth.blogspot.dk/2014/10/easy-clustered-standard-errors-in-r.html

De la documentation:

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)

Chandler
la source

pour moi coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0")) ainsi que coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0",cluster="group"))produire des résultats identiques. Savez-vous pourquoi c'est le cas?

Peter Pan

1

Le lien people.su.se/~ma/clustering.pdf ne fonctionne plus. Vous souvenez-vous du titre de la page?

MERose

8

Après beaucoup de lecture, j'ai trouvé la solution pour la mise en cluster dans le lmcadre.

Mahmood Arai a publié un excellent livre blanc qui fournit un didacticiel sur la mise en cluster dans le lmcadre, qu'il fait avec des corrections de degrés de liberté au lieu de mes tentatives en désordre décrites ci-dessus. Il fournit ses fonctions pour les deux à une et deux voies matrices de covariance de regroupement ici .

Enfin, bien que le contenu ne soit pas disponible gratuitement, le document Principalement inoffensif, Econometrics d’ Angrist et Pischke, contient une section très utile sur la mise en grappes.

Mise à jour du 27/04/2015 pour ajouter le code d'un article de 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
la source

1

Le papier d'Arai n'est plus en ligne. Pouvez-vous fournir le lien réel?

MERose

@MERose - Désolé! Malheureusement, nous ne pouvons pas joindre des fichiers PDF! J'ai trouvé cet article de blog qui compare le code. Je vais éditer cette réponse pour inclure le code.

Richard Herron

Cela pourrait être une version mise à jour du livre blanc: Arai Mahmood, erreurs standards Cluster robustes à l' aide de R .

gung - Réintégrer Monica

4

Le moyen le plus simple de calculer les erreurs standard de cluster dans R consiste à utiliser la fonction récapitulative modifiée.

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

Il y a un excellent article sur le clustering dans le cadre de lm. Le site fournit également la fonction de résumé modifiée pour les clusters à une et à deux voies. Vous pouvez trouver la fonction et le tutoriel ici .

Tolga Suer
la source

1

Si vous n'avez pas d' timeindex, vous n'en avez pas besoin: vous plmen ajouterez un fictif, qui ne sera utilisé que si vous le demandez. Donc, cet appel devrait fonctionner :

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

Sauf que ce n'est pas le cas, ce qui suggère que vous avez rencontré un bogue plm. (Ce bogue a maintenant été corrigé dans SVN. Vous pouvez installer la version de développement à partir d'ici .)

Mais puisque ce serait de toute façon un timeindex fictif , nous pouvons le créer nous-mêmes:

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

Maintenant cela fonctionne:

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

Remarque importante : vcovHC.plm()in plmestimera par défaut Arellano en cluster par groupe de SE. Ce qui diffère de l'estimation vcovHC.lm()in sandwich(par exemple, les vcovHCSE de la question initiale), c'est-à-dire des SE à cohérence d'hétéroscédasticité sans regroupement.

Une approche distincte consiste à s'en tenir aux lmrégressions de variables nominales et au package 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

Dans les deux cas, vous obtiendrez les SE d'Arellano (1987) avec regroupement par groupe. Ce multiwayvcovpaquet est une évolution directe et significative des fonctions de clustering originales d’Arai.

Vous pouvez également consulter la matrice de variance-covariance résultante des deux approches, donnant la même estimation de variance pour 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
la source

Veuillez consulter ce lien: multiwayvcov est amorti: sites.google.com/site/npgraham1/research/code

HoneyBuddha

Regroupement des erreurs standard dans R (manuellement ou dans plm)

Réponses: