R中负二项回归的稳健标准误差与Stata的标准误差不匹配

Question

我正在R中复制负二项回归模型。在计算稳健标准误差时，输出与标准误差的Stata输出不匹配。

Stata的原始代码是

nbreg displaced  eei lcostofwar cfughh roadskm lpopdensity ltkilled, robust nolog

我尝试了手动计算和sandwich中的vcovHC。然而，两者都没有产生相同的结果。

我的回归模型如下：mod1 <- glm.nb(displaced ~ eei + costofwar_log + cfughh + roadskm + popdensity_log + tkilled_log, data = mod1_df)

使用vcovHC，我尝试了从HC0到HC5的所有选项。尝试1:

cov_m1 <- vcovHC(mod1, type = "HC0", sandwich = T)
se <- sqrt(diag(cov_m1))

尝试2:

mod1_rob <- coeftest(mod1, vcovHC = vcov(mod1, type = "HC0"))

最成功的是HC0和vcov = sandwich，但没有任何SE是正确的。

有什么建议吗？

编辑

我的输出如下(使用HC0(：

Estimate Std. Error z value  Pr(>|z|)    
(Intercept)     1.3281183  1.5441312  0.8601  0.389730    
eei            -0.0435529  0.0183359 -2.3753  0.017536 *  
costofwar_log   0.2984376  0.1350518  2.2098  0.027119 *  
cfughh         -0.0380690  0.0130254 -2.9227  0.003470 ** 
roadskm         0.0020812  0.0010864  1.9156  0.055421 .  
popdensity_log -0.4661079  0.1748682 -2.6655  0.007688 ** 
tkilled_log     1.0949084  0.2159161  5.0710 3.958e-07 ***

我试图复制的Stata输出是：

Estimate Std. Error    
(Intercept)     1.328     1.272
eei            -0.044     0.015 
costofwar_log   0.298     0.123  
cfughh         -0.038     0.018 
roadskm         0.002     0.0001   
popdensity_log -0.466     0.208 
tkilled_log     1.095     0.209  

数据集位于此处，记录的变量为：

mod1_df <- table %>% 
select(displaced, eei_01, costofwar, cfughh, roadskm, popdensity, 
tkilled)
mod1_df$popdensity_log <- log(mod1_df$popdensity + 1)
mod1_df$tkilled_log <- log(mod1_df$tkilled + 1)
mod1_df$costofwar_log <- log(mod1_df$costofwar + 1)
mod1_df$eei <- mod1_df$eei_01*100

Answer 1

Stata使用观测到的Hessian进行计算，glm.nb()使用预期的Hessin。因此，sandwich()函数使用的默认bread()是不同的，从而导致不同的结果。还有其他R包使用观测到的hessian来进行其方差-协方差估计(例如，gamlss(，但这些包没有为sandwich包提供estfun()方法。

因此，下面我简单地设置了一个专用的bread_obs()函数，该函数从negbin对象中提取ML估计，设置负对数似然，通过numDeriv::hessian()数值计算观测到的Hessian，并计算"；面包"；从中(省略对数(θ(的估计(：

bread_obs <- function(object, method = "BFGS", maxit = 5000, reltol = 1e-12, ...) {
## data and estimated parameters
Y <- model.response(model.frame(object))
X <- model.matrix(object)
par <- c(coef(object), "log(theta)" = log(object$theta))
## dimensions
n <- NROW(X)
k <- length(par)
## nb log-likelihood
nll <- function(par) suppressWarnings(-sum(dnbinom(Y,
mu = as.vector(exp(X %*% head(par, -1))),
size = exp(tail(par, 1)), log = TRUE)))
## covariance based on observed Hessian
rval <- numDeriv::hessian(nll, par) 
rval <- solve(rval) * n
rval[-k, -k]
}

使用该函数，我可以将sandwich()输出(基于预期的Hessian(与使用bread_obs()的输出(基于观察到的Hessin(进行比较。

s_exp <- sandwich(mod1)
s_obs <- sandwich(mod1, bread = bread_obs)
cbind("Coef" = coef(mod1), "SE (Exp)" = sqrt(diag(s_exp)), "SE (Obs)" = sqrt(diag(s_obs)))
##                  Coef SE (Exp) SE (Obs)
## (Intercept)     1.328    1.259    1.259
## eei            -0.044    0.017    0.015
## costofwar_log   0.298    0.160    0.121
## cfughh         -0.038    0.015    0.018
## roadskm         0.002    0.001    0.001
## popdensity_log -0.466    0.135    0.207
## tkilled_log     1.095    0.179    0.208

与Stata相比，这仍然有轻微的差异，但这些可能是与优化等的数值差异。

如果为negbin对象创建新的专用bread()方法

bread.negbin <- bread_obs

如果执行CCD_ 23，则方法dispatch将使用此方法。

Answer 2

在R中，您需要手动提供自由度校正，所以请尝试我从以下来源借用的方法：

dfa <- (G/(G - 1)) * (N - 1)/pm1$df.residual
# display with cluster VCE and df-adjustment
firm_c_vcov <- dfa * vcovHC(pm1, type = "HC0", cluster = "group", adjust = T)
coeftest(pm1, vcov = firm_c_vcov)

这里，G是数据集中面板的数量，N是观测的数量，而pm1是模型估计的数量。显然，您可以放弃集群。