我正在 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
最佳答案
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 将使用它。
关于R 中负二项式回归的稳健标准误差与 Stata 中的不匹配,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/53629037/