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
). 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!
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
—
mpiktas,
clusteropcja jest wyjaśniona. Jestem pewien, że replikacja w R. byłaby możliwa
+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,
—
hans0l0
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.
@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