0

我正在尝试重现此stata示例并从. 数据可在此处获得。stargazertexreg

要运行回归并获得 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 值？

4

2 回答 2

4

强大的标准错误texreg很容易：只需直接通过 coeftest！

自从上次回答问题以来，这变得容易多了：看来您现在可以直接通过带有所需方差 - 协方差矩阵的 coeftest。缺点：您会失去拟合统计的优度（例如 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)
于 2018-02-27T18:40:07.520 回答
2

首先，请注意，as.integer一旦使用具有非数字行名的数据，您的使用是危险的并且可能会导致问题。例如，使用行名由汽车名称组成的内置数据集mtcars，您的函数会将所有行名强制转换为NA，并且您的函数将不起作用。

对于您的实际问题，您可以为 提供自定义 p 值texreg，这意味着您需要计算相应的 p 值。为此，您可以计算方差-协方差矩阵，计算检验统计量，然后手动计算 p 值，或者您只需计算方差-协方差矩阵并将其提供给 eg 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 分布的内置函数（参见 eg ?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
于 2016-03-29T14:20:51.927 回答