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

Question

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

原始的Stata代码是

nbreg displaced  eei lcostofwar cfughh roadskm lpopdensity ltkilled, robust nolog

我尝试了手动计算和 vcovHC来自 sandwich .但是，两者都不会产生相同的结果。

我的回归模型如下: 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  

找到数据集here重新编码的变量是:

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() 使用预期的 Hessian 矩阵。因此，sandwich()函数使用的默认bread()是不同的，导致不同的结果。还有其他 R 包使用观察到的 hessian 来估计其方差-协方差(例如，gamlss)，但这些包没有为 提供 estfun() 方法三明治包。

因此，下面我简单地设置了一个专用的 bread_obs() 函数，它从 negbin 对象中提取 ML 估计值，设置负对数似然，计算通过 numDeriv::hessian() 以数字方式观察 Hessian 并从中计算“面包”(省略对 log(theta) 的估计):

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() 的输出(基于观察到的 Hessian)进行比较.

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

如果你执行 sandwich(mod1)，方法 dispatch 将使用它。