texreg 的集群标准错误？
Clustered standard errors with texreg?
我正在尝试重现 this stata example and move from stargazer to texreg. The data is available here。
到运行回归并得到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 值？
首先，请注意您对 as.integer 的使用是危险的，并且一旦您使用具有非数字行名的数据可能会导致问题。例如，使用行名由汽车名称组成的内置数据集 mtcars，您的函数会将所有行名强制转换为 NA，并且您的函数将不起作用。
对于您的实际问题，您可以向texreg提供自定义p值，这意味着您需要计算相应的p值。为此，您可以计算方差-协方差矩阵，计算检验统计量，然后手动计算 p 值，或者您只需计算方差-协方差矩阵并将其提供给例如coeftest。然后您可以从那里提取标准误差和 p 值。由于我不愿意下载任何数据，因此我将 mtcars-data 用于以下内容：
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 的稳健标准错误很简单：直接通过系数测试即可！
自上次回答问题以来，这变得容易多了：看来您现在可以直接通过所需 variance-covariance 矩阵的系数测试。缺点：你失去了拟合优度统计（例如 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)
我正在尝试重现 this stata example and move from stargazer to texreg. The data is available here。
到运行回归并得到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 值？
首先，请注意您对 as.integer 的使用是危险的，并且一旦您使用具有非数字行名的数据可能会导致问题。例如，使用行名由汽车名称组成的内置数据集 mtcars，您的函数会将所有行名强制转换为 NA，并且您的函数将不起作用。
对于您的实际问题，您可以向texreg提供自定义p值，这意味着您需要计算相应的p值。为此，您可以计算方差-协方差矩阵，计算检验统计量，然后手动计算 p 值，或者您只需计算方差-协方差矩阵并将其提供给例如coeftest。然后您可以从那里提取标准误差和 p 值。由于我不愿意下载任何数据，因此我将 mtcars-data 用于以下内容：
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 的稳健标准错误很简单：直接通过系数测试即可！
自上次回答问题以来，这变得容易多了：看来您现在可以直接通过所需 variance-covariance 矩阵的系数测试。缺点：你失去了拟合优度统计（例如 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)