如何解决如何并排比较零膨胀负二项式回归输出，有和没有聚集错误
在 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),dist = "negbin")
v_lake = vcovCL(with.clusters,type = "HC1",cluster = ~lake)
with.clusters.final <- coeftest(with.clusters,v_lake)
如何生成比较对象 no.cluster 和 with.clusters.final 的表格？
解决方法
我意识到我错误地考虑了这个问题。与其直接将 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),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,"t")),+ p = list(NULL,"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
这里有一个更简单的方法：对无簇对象使用 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 参数来选择行并为其命名。
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),dist = "negbin")
modelsummary(list(mod1,mod2),vcov = list("classical",~lake),stars = TRUE)
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
