R中具有稳健聚集标准误差的逻辑回归答案

【问题标题】：Logistic regression with robust clustered standard errors in RR中具有稳健聚集标准误差的逻辑回归
【发布时间】：2013-05-06 02:48:58
【问题描述】：

一个新手问题：有人知道如何在 R 中使用聚集标准错误运行逻辑回归吗？在 Stata 中它只是 logit Y X1 X2 X3, vce(cluster Z)，但不幸的是我还没有弄清楚如何在 R 中进行相同的分析。提前致谢！

【问题讨论】：

vcovHC() 包中的 vcovHC() 函数也可能有用（不确定它是否适用于逻辑回归估计）
如果您从 Stata 迁移，您可能会发现名为 plm 的软件包很有用。此外，还有一个名为 pcse 的包，用于通过在估计后操纵方差协方差矩阵来实现面板校正标准误差
非常感谢您的回复，本和休伯特。我还将测试您建议的软件包，看看它们是否适用于逻辑估计。再次感谢！

标签： r regression stata

【解决方案1】：

您可能想查看rms（回归建模策略）包。所以，lrm 是逻辑回归模型，如果 fit 是你的输出名称，你会有这样的东西：

fit=lrm(disease ~ age + study + rcs(bmi,3), x=T, y=T, data=dataf)

fit

robcov(fit, cluster=dataf$id)

bootcov(fit,cluster=dataf$id)

您必须在模型语句中指定x=T、y=T。 rcs 表示 3 节的受限三次样条。

【讨论】：

非常感谢！它创造了奇迹！我会更仔细地阅读 rms 的手册，看看是否有办法按国家和年份对系数进行聚类。再次感谢您！
这个答案已经非常好，但如果它完全可复制，则可以改进。我不知道变量来自哪里，输出是什么，以及为什么需要rcs(bmi,3)。

【解决方案2】：

过去两天我一直在努力解决这个问题；我神奇地发现了一个似乎注定要实现伟大事物的新包——例如，我还在分析中运行了一些集群健壮的 Tobit 模型，并且这个包也内置了该功能。更不用说语法比我见过的所有其他解决方案都更干净（我们说的是接近 Stata 级别的干净）。

所以对于你的玩具示例，我会运行：

library(Zelig)
logit<-zelig(Y~X1+X2+X3,data=data,model="logit",robust=T,cluster="Z")

等等！

【讨论】：

哇，这似乎确实以我的 R 代码从未出现过的方式“正常工作”。这是新功能吗？如果不是，为什么 Zelig 不是在 R 中解决这个问题的规范方法？
不知道，但我希望是这样。 The project 看起来确实雄心勃勃！ [Google 群组](groups.google.com/forum/m/#!forum/zelig-statistical-software) 似乎不太活跃，所以不确定进展有多快。
不幸的是，我认为该命令在最新版本的Zelig（在 CRAN 上）中不起作用。我刚刚运行了几个带有和不带有 cluster 参数的模型，标准错误完全相同。我相信从 4.0 版开始就是这样，我最后一次使用这个包。
是的，他们确实暂时放弃了该功能。检查他们的谷歌小组（去他们网站的社区部分）——他们正在重组整个项目；一位开发人员在回复我的帖子时说，他们正在努力恢复集群/强大的功能
大约三年后，集群功能未恢复：Error in glm.control(cluster = "group") : unused argument (cluster = "group")。

【解决方案3】：

R 包miceadds 中有一个命令glm.cluster，它似乎给出了与Stata 使用选项vce(cluster) 相同的逻辑回归结果。请参阅文档here。

在本页的一个示例中，命令

mod2 <- miceadds::glm.cluster(data=dat, formula=highmath ~ hisei + female,
                              cluster="idschool", family="binomial")
summary(mod2)

提供与 Stata 命令相同的稳健标准错误

logit highmath hisei female, vce(cluster idschool)

例如变量 hisei 的标准误差为 0.004038。

【讨论】：

【解决方案4】：

另一种选择是使用sandwich 和lmtest 包，如下所示。假设z 是数据集中dat 中包含聚类指标的列。那么

# load libraries
library("sandwich")
library("lmtest")

# fit the logistic regression
fit = glm(y ~ x, data = dat, family = binomial)

# get results with clustered standard errors (of type HC0)
coeftest(fit, vcov. = vcovCL(fit, cluster = dat$z, type = "HC0"))

会做的。

【讨论】：