如何并排比较零膨胀负二项式回归输出，有和没有聚集错误答案

【问题标题】：how to compare zero-inflated negative binomial regression outputs side-by-side, with and without clustered errors如何并排比较零膨胀负二项式回归输出，有和没有聚集错误
【发布时间】：2021-08-21 21:00:19
【问题描述】：

在 R 的并排表中比较回归模型的方法有很多，包括包stargazer、huxtable 和gtsummary。我正在努力使用两个零膨胀负二项式模型来做到这一点，其中一个具有聚集标准误差而另一个没有。

问题在于不同的对象类别。聚集错误模型属于pscl 包中的“coeftest”类，而非聚集模型属于lmtest 包中的“zeroinfl”类。

这是一个简单的例子。

library(tidyverse)
library(pscl)
library(sandwich)
library(lmtest)

## data
data("bioChemists", package = "pscl")

set.seed(3.14)
zinb <- read.csv("https://stats.idre.ucla.edu/stat/data/fish.csv") %>%
  # add some random data representing lakes
  mutate(lake = sample(c("lake1", "lake2", "lake3"), n(), replace = TRUE))

# zero inflated model without clusters
no.clusters <- zeroinfl(count ~ child + camper,
                        data = zinb, dist = "negbin")

# cluster by lake
with.clusters <- zeroinfl(count ~ child + camper + factor(lake),
                          data = zinb, dist = "negbin")
v_lake = vcovCL(with.clusters, type = "HC1", cluster = ~lake)
with.clusters.final <- coeftest(with.clusters, v_lake)

如何生成一个比较对象no.cluster 和with.clusters.final 的表格？

【问题讨论】：

我从vcovCL 电话中收到警告消息“警告消息：在 Ops.factor(child + camper, factor(lake)) 中：'+' 对因素没有意义”。你明白为什么吗？（也许与您的问题无关。）
@BenBolker，我对此进行了调查，只有在指定第二个模型公式（零通胀 logit 模型）时才会出现该警告消息。我不明白这种行为，但我已编辑代码以删除警告。

标签： r stargazer

【解决方案1】：

我意识到我在错误地思考这个问题。我无需尝试直接将no.cluster 与使用coeftest 生成的对象进行比较，而是只需将with.clusters 中的标准错误替换为coeftest 生成的稳健值。

这比听起来要复杂一些，因为我使用的是零膨胀回归。

零膨胀回归由两个模型组成，这两个模型都包含在coeftest 的输出中。我只关心以count_开头的值。

> with.clusters.final

t test of coefficients:

                         Estimate Std. Error t value              Pr(>|t|)    
count_(Intercept)         1.08748    0.11608  9.3686 < 0.00000000000000022 ***
count_child              -0.82116    0.33572 -2.4460             0.0151678 *  
count_camper              0.80562    0.15933  5.0563        0.000000850945 ***
count_factor(lake)lake2   0.18206    0.19467  0.9352             0.3506175    
count_factor(lake)lake3  -0.59536    0.10789 -5.5184        0.000000088810 ***
zero_(Intercept)        -21.32602    5.36229 -3.9770        0.000092516485 ***
zero_child               11.73360    2.42002  4.8486        0.000002241063 ***
zero_camper              -3.07137    0.85151 -3.6070             0.0003772 ***
zero_factor(lake)lake2   10.87828    2.60947  4.1688        0.000042846049 ***
zero_factor(lake)lake3    1.35318    0.22321  6.0624        0.000000005192 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我需要删除以zero_ 开头的行，然后从其余名称中删除“count_”前缀。

> clustered.se <- with.clusters.final[,2]
> clustered.se <- clustered.se[str_sub(names(clustered.se), 1, 5) == "count"]
> names(clustered.se) <- str_remove(names(clustered.se), "count_")
> clustered.se
      (Intercept)             child            camper factor(lake)lake2 factor(lake)lake3 
        0.1160763         0.3357177         0.1593300         0.1946651         0.1078864

我们可以对检验统计量和 p 值做同样的事情。这是一个简化流程的函数。

clustered_stat <- function(coeftest_object, statistic){
  statistic.row <- if(statistic == "se"){
    2
  } else if(statistic == "t") {
    3
  } else if(statistic == "p") {
    4
  }
  
  # pull the appropraite row out of the coeftest object
  new.values <- coeftest_object[,statistic.row]
  # just keep the values from the count model, not the logit model
  new.values <- new.values[str_sub(names(new.values), 1, 5) == "count"]
  # remove the "count_" prefix from the names of the values
  names(new.values) <- str_remove(names(new.values), "count_")
  new.values
}

现在我可以将原始的集群 zeroinfl 对象提供给 stargazer，用 coeftest 创建的强大版本替换统计信息。

> stargazer::stargazer(no.clusters, with.clusters,
+                      se = list(NULL, clustered_stat(with.clusters.final, "se")),
+                      t = list(NULL, clustered_stat(with.clusters.final, "t")),
+                      p = list(NULL, clustered_stat(with.clusters.final, "p")),
+                      type = "text",
+                      column.labels = c("no clusters", "with clusters"))

==============================================
                      Dependent variable:     
                  ----------------------------
                             count            
                   no clusters   with clusters
                       (1)            (2)     
----------------------------------------------
child               -0.911***      -0.821**   
                     (0.285)        (0.336)   
                                              
camper               0.798***      0.806***   
                     (0.305)        (0.159)   
                                              
factor(lake)lake2                    0.182    
                                    (0.195)   
                                              
factor(lake)lake3                  -0.595***  
                                    (0.108)   
                                              
Constant             1.052***      1.087***   
                     (0.270)        (0.116)   
                                              
----------------------------------------------
Observations           250            250     
Log Likelihood       -434.906      -430.784   
==============================================
Note:              *p<0.1; **p<0.05; ***p<0.01

【讨论】：

下面更简单的方法。

【解决方案2】：

这里有一个更简单的方法：将coeftest 用于无集群对象，但没有集群。然后你可以使用例如huxtable::huxreg() 打印出你的系数：

with.clusters.final <- coeftest(with.clusters, v_lake, save = TRUE)
no.clusters.coeftest <- coeftest(no.clusters, save = TRUE)
huxreg(no.clusters.coeftest, with.clusters.final)
                        ───────────────────────────────────────────────────────
                                                            (1)           (2)  
                                                  ─────────────────────────────
                          count_(Intercept)           1.052 ***     1.087 ***  
                                                     (0.270)       (0.116)     
                          count_child                -0.911 **     -0.821 *    
                                                     (0.285)       (0.336)     
                          count_camper                0.798 **      0.806 ***  
                                                     (0.305)       (0.159)     
                          zero_(Intercept)          -11.499       -21.326 ***  
                                                    (55.699)       (5.362)     
                          zero_child                 10.483        11.734 ***  
                                                    (55.659)       (2.420)     
                          zero_camper                -9.501        -3.071 ***  
                                                    (55.663)       (0.852)     
                          count_factor(lake)lake2                   0.182      
                                                                   (0.195)     
                          count_factor(lake)lake3                  -0.595 ***  
                                                                   (0.108)     
                          zero_factor(lake)lake2                   10.878 ***  
                                                                   (2.609)     
                          zero_factor(lake)lake3                    1.353 ***  
                                                                   (0.223)     
                                                  ─────────────────────────────
                          nobs                                                 
                        ───────────────────────────────────────────────────────
                          *** p < 0.001; ** p < 0.01; * p < 0.05.

可以使用coefs 参数来选择行并为其命名。

【讨论】：

您可能还想告诉 Ben Bolker 将 zeroinfl 和 hurdle 模型包含到 broom.mixed....

【解决方案3】：

The modelsummary package 接受 vcov 参数，该参数可以是字符串列表（例如，"robust" 或 "classical"）或带有集群变量的单边公式。在后台，它会自动调用 sandwich 包中的函数来计算不确定性估计。（免责声明：我是包的维护者）。

请注意，表格底部会自动附加一个注释以指示标准错误的种类。

library(pscl)
library(modelsummary)

set.seed(3.14)
zinb <- read.csv("https://stats.idre.ucla.edu/stat/data/fish.csv")
zinb$lake <- sample(c("lake1", "lake2", "lake3"), nrow(zinb), replace = TRUE)

mod1 <- zeroinfl(count ~ child + camper,
                 data = zinb, dist = "negbin")

mod2 <- zeroinfl(count ~ child + camper + factor(lake),
                 data = zinb, dist = "negbin")

modelsummary(list(mod1, mod2),
             vcov = list("classical", ~lake),
             stars = TRUE)

【讨论】：