2

我在 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

数据集位于此处，重新编码的变量为：

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
4

2 回答 2

3

Stata 使用观察到的 Hessian 进行计算，glm.nb()使用预期的 Hessian。bread()因此，函数采用的默认值sandwich()不同，导致结果不同。还有其他 R 包使用观察到的 hessian 进行方差协方差估计（例如，gamlss），但这些包不提供包的estfun()方法sandwich。

因此，下面我简单地设置了一个专用bread_obs()函数，该函数从对象中提取 ML 估计negbin，设置负对数似然，通过数值计算观察到的 HessiannumDeriv::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, vcov = 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 相比，这仍然存在细微差异，但这些可能是优化等方面的数值差异。

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

bread.negbin <- bread_obs

那么如果你这样做，方法调度将使用它sandwich(mod1)。

于 2019-01-16T02:18:56.517 回答
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是您的模型估计的。显然，您可以删除集群。

于 2018-12-05T10:25:14.483 回答