Кластерные стандартные ошибки с texreg?
Я пытаюсь воспроизвести этот пример и перейти от stargazer в texreg, Данные доступны здесь.
Чтобы запустить регрессию и получить SE я запускаю этот код:
library(readstata13)
library(sandwich)
cluster_se <- function(model_result, data, cluster){
model_variables <- intersect(colnames(data), c(colnames(model_result$model), cluster))
model_rows <- as.integer(rownames(model_result$model))
data <- data[model_rows, model_variables]
cl <- data[[cluster]]
M <- length(unique(cl))
N <- nrow(data)
K <- model_result$rank
dfc <- (M/(M-1))*((N-1)/(N-K))
uj <- apply(estfun(model_result), 2, function(x) tapply(x, cl, sum));
vcovCL <- dfc*sandwich(model_result, meat=crossprod(uj)/N)
sqrt(diag(vcovCL))
}
elemapi2 <- read.dta13(file = 'elemapi2.dta')
lm1 <- lm(formula = api00 ~ acs_k3 + acs_46 + full + enroll, data = elemapi2)
se.lm1 <- cluster_se(model_result = lm1, data = elemapi2, cluster = "dnum")
stargazer::stargazer(lm1, type = "text", style = "aer", se = list(se.lm1))
==========================================================
api00
----------------------------------------------------------
acs_k3 6.954
(6.901)
acs_46 5.966**
(2.531)
full 4.668***
(0.703)
enroll -0.106**
(0.043)
Constant -5.200
(121.786)
Observations 395
R2 0.385
Adjusted R2 0.379
Residual Std. Error 112.198 (df = 390)
F Statistic 61.006*** (df = 4; 390)
----------------------------------------------------------
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.
texreg производит это:
texreg::screenreg(lm1, override.se=list(se.lm1))
========================
Model 1
------------------------
(Intercept) -5.20
(121.79)
acs_k3 6.95
(6.90)
acs_46 5.97 ***
(2.53)
full 4.67 ***
(0.70)
enroll -0.11 ***
(0.04)
------------------------
R^2 0.38
Adj. R^2 0.38
Num. obs. 395
RMSE 112.20
========================
Как я могу исправить p-значения?
2 ответа
Во-первых, обратите внимание, что ваше использование as.integer опасно и может вызвать проблемы при использовании данных с нечисловыми именами строк. Например, используя встроенный набор данных mtcars чьи имена строк состоят из имен машин, ваша функция приведет все имена строк к NA, и ваша функция не будет работать.
На ваш актуальный вопрос, вы можете предоставить пользовательские p-значения для texreg, что означает, что вам нужно вычислить соответствующие p-значения. Чтобы достичь этого, вы можете вычислить матрицу дисперсии-ковариации, вычислить статистику теста, а затем вручную вычислить значение p или просто вычислить матрицу дисперсии-ковариации и предоставить ее, например, coeftest, Затем вы можете извлечь стандартные ошибки и p-значения оттуда. Поскольку я не желаю загружать какие-либо данные, я использую mtcars-данные для следующего:
library(sandwich)
library(lmtest)
library(texreg)
cluster_se <- function(model_result, data, cluster){
model_variables <- intersect(colnames(data), c(colnames(model_result$model), cluster))
model_rows <- rownames(model_result$model) # changed to be able to work with mtcars, not tested with other data
data <- data[model_rows, model_variables]
cl <- data[[cluster]]
M <- length(unique(cl))
N <- nrow(data)
K <- model_result$rank
dfc <- (M/(M-1))*((N-1)/(N-K))
uj <- apply(estfun(model_result), 2, function(x) tapply(x, cl, sum));
vcovCL <- dfc*sandwich(model_result, meat=crossprod(uj)/N)
}
lm1 <- lm(formula = mpg ~ cyl + disp, data = mtcars)
vcov.lm1 <- cluster_se(model_result = lm1, data = mtcars, cluster = "carb")
standard.errors <- coeftest(lm1, vcov. = vcov.lm1)[,2]
p.values <- coeftest(lm1, vcov. = vcov.lm1)[,4]
texreg::screenreg(lm1, override.se=standard.errors, override.p = p.values)
И для полноты картины давайте сделаем это вручную:
t.stats <- abs(coefficients(lm1) / sqrt(diag(vcov.lm1)))
t.stats
(Intercept) cyl disp
38.681699 5.365107 3.745143
Это ваша t-статистика с использованием устойчивых к кластеру стандартных ошибок. Степень свободы хранится в lm1$df.residualи использование встроенных функций для t-распределения (см., например, ?pt), мы получаем:
manual.p <- 2*pt(-t.stats, df=lm1$df.residual)
manual.p
(Intercept) cyl disp
1.648628e-26 9.197470e-06 7.954759e-04
Вот, pt является функцией распределения, и мы хотим вычислить вероятность наблюдения статистики, по крайней мере, такой же экстремальной, как та, которую мы наблюдаем. Поскольку мы тестируем двустороннюю и это симметричную плотность, мы сначала берем левую крайность, используя отрицательное значение, а затем удваиваем его. Это идентично использованию 2*(1-pt(t.stats, df=lm1$df.residual)), Теперь просто проверьте, что это дает тот же результат, что и раньше:
all.equal(p.values, manual.p)
[1] TRUE
Надежные стандартные ошибки с texreg это просто: просто пройдите самый прямой!
Это стало намного проще с тех пор, как на вопрос был дан последний ответ: кажется, что теперь вы можете просто пройти самое сложное с нужной матрицей дисперсии-ковариации напрямую. Недостаток: вы теряете достоверность статистики соответствия (например, R^2 и количество наблюдений), но в зависимости от ваших потребностей это может не быть большой проблемой
Как включить надежные стандартные ошибки с texreg
> screenreg(list(reg1, coeftest(reg1,vcov = vcovHC(reg1, 'HC1'))),
custom.model.names = c('Standard Standard Errors', 'Robust Standard Errors'))
=============================================================
Standard Standard Errors Robust Standard Errors
-------------------------------------------------------------
(Intercept) -192.89 *** -192.89 *
(55.59) (75.38)
x 2.84 ** 2.84 **
(0.96) (1.04)
-------------------------------------------------------------
R^2 0.08
Adj. R^2 0.07
Num. obs. 100
RMSE 275.88
=============================================================
*** p < 0.001, ** p < 0.01, * p < 0.05
Чтобы сгенерировать этот пример, я создал фрейм данных с гетероскедастичностью, полный пример кода приведен ниже:
require(sandwich);
require(texreg);
set.seed(1234)
df <- data.frame(x = 1:100);
df$y <- 1 + 0.5*df$x + 5*100:1*rnorm(100)
reg1 <- lm(y ~ x, data = df)