Надежные стандартные ошибки для отрицательной биномиальной регрессии в 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  

Набор данных найден здесь , а перекодированные переменные:

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 использует наблюдаемый гессиан для своих вычислений, glm.nb() использует ожидаемый гессиан.Поэтому значение по умолчанию bread(), используемое функцией sandwich(), отличается, что приводит к разным результатам.Существуют другие пакеты R, которые используют наблюдаемый гессиан для оценки дисперсии-ковариации (например, gamlss), но они не предоставляют метод estfun() для пакета sandwich.

Следовательно, ниже Iпросто установите выделенную функцию bread_obs(), которая извлекает оценки ML из объекта negbin, устанавливает отрицательное логарифмическое правдоподобие, численно вычисляет наблюдаемый гессиан через numDeriv::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() (на основе ожидаемого гессиана) с выводом, используя bread_obs() (на основе наблюдаемого гессиана).

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).

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 - ваша оценочная модель.Очевидно, вы можете отказаться от кластеризации.