在R中进行Probit和Logit回归的稳健聚类标准误差

Question

在R中进行Probit和Logit回归的稳健聚类标准误差

8

最初，我主要想在R中使用带聚类标准误的probit / logit模型，这在Stata中非常直观。我在这里找到了答案 Logistic regression with robust clustered standard errors in R。

因此，我尝试比较使用鲁棒标准误和聚类标准误的Stata和R的结果。但是我注意到两种软件的标准误输出并不完全相同。然而，如果我使用这里建议的方法https://diffuseprior.wordpress.com/2012/06/15/standard-robust-and-clustered-standard-errors-computed-in-r/。我可以从R和Stata获得线性回归的精确输出。因此，我担心我在R中编写的代码是否正确，以及如果我想运行probit模型而不是logit模型应该使用什么命令。或者是否有任何优雅的替代方案来解决这个问题？谢谢。

R代码

## 1. linear regression
library(rms) 
# model<-lm(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width,iris)
summary(model)
fit=ols(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width, x=T, y=T, data=iris)
fit
robcov(fit) #robust standard error
robcov(fit, cluster=iris$Species) #clustered standard error


## 2. logistic regression
##demo data generation   
set.seed(1234)
subj<-rep(1:20,each=4)
con1<-rep(c(1,0),40)
con2<-rep(c(1,1,0,0),20) 
effect<-rbinom(80,1,0.34)
data<-data.frame(subj,con1,con2,effect)
library(foreign);write.dta(data,'demo_data.dta')

library(rms)
fit=lrm(effect ~ con1 + con2, x=T, y=T, data=data)
fit
robcov(fit)  ##robust standard error
robcov(fit, cluster=data$subj) ## clustered standard error

Stata 代码

## 1. linear regression
webuse iris
reg seplen sepwid petlen petwid
reg seplen sepwid petlen petwid,r
reg seplen sepwid petlen petwid,cluster(iris)


## 2. logistic regression

use demo_data,clear
logit effect con1 con2
logit effect con1 con2,r
logit effect con1 con2,cluster(subj)

- johnsonzhj

2

你能具体说明一下 not exactly the same 是什么意思吗？这里涉及到很多默认值，它们可能是不同的。我们无法事先确定哪些默认值更好。但如果你想得到完全相同的值，你需要弄清楚 Stata 和 robcov 使用了哪些默认值，并进行相应的调整。 - coffeinjunky

感谢您的评论，我已经编辑了我的问题以提供更多信息。 - johnsonzhj

你是否可能在运行 logistic 前就使用了 logit？“logistic 显示估计值为比率；若要查看系数，请在运行 logistic 后输入 logit”。 - noumenal

Stata的输出与通过vcovHC获得的输出相同。实际上，我刚刚检查了这里提出的方法：链接，我可以获得与从stata获得的输出相同的输出，所以仍然很好奇如何使用rms包来实现这一点。 - johnsonzhj

知道了，谢谢，我会仔细阅读并稍后给你反馈。 - johnsonzhj

显示剩余3条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- coffeinjunky · Accepted Answer

我更喜欢使用sandwich包来计算鲁棒标准误。其中一个原因是其出色的文档。请参见vignette("sandwich")，其中清楚地显示了所有可用的默认值和选项，以及相应文章，该文章解释了如何使用自定义bread和meat处理特殊情况下的?sandwich。

我们可以使用sandwich来找出您发布的选项之间的差异。差异很可能是自由度校正。以下是简单线性回归的比较：

library(rms)
library(sandwich)

fitlm <-lm(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width,iris)

#Your Blog Post:
X <- model.matrix(fitlm)
n <- dim(X)[1]; k <- dim(X)[2]; dfc <- n/(n-k)    
u <- matrix(resid(fitlm))
meat1 <- t(X) %*% diag(diag(crossprod(t(u)))) %*% X
Blog <- sqrt(dfc*diag(solve(crossprod(X)) %*% meat1 %*% solve(crossprod(X))))

# rms fits:
fitols <- ols(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width, x=T, y=T, data=iris)
Harrell <- sqrt(diag(robcov(fitols, method="huber")$var))
Harrell_2 <- sqrt(diag(robcov(fitols, method="efron")$var))

# Variations available in sandwich:    
variations <- c("const", "HC0", "HC1", "HC2","HC3", "HC4", "HC4m", "HC5")
Zeileis <- t(sapply(variations, function(x) sqrt(diag(vcovHC(fitlm, type = x)))))
rbind(Zeileis, Harrell, Harrell_2, Blog)

          (Intercept) Sepal.Width Petal.Length Petal.Width
const       0.2507771  0.06664739   0.05671929   0.1275479
HC0         0.2228915  0.05965267   0.06134461   0.1421440
HC1         0.2259241  0.06046431   0.06217926   0.1440781
HC2         0.2275785  0.06087143   0.06277905   0.1454783
HC3         0.2324199  0.06212735   0.06426019   0.1489170
HC4         0.2323253  0.06196108   0.06430852   0.1488708
HC4m        0.2339698  0.06253635   0.06482791   0.1502751
HC5         0.2274557  0.06077326   0.06279005   0.1454329
Harrell     0.2228915  0.05965267   0.06134461   0.1421440
Harrell_2   0.2324199  0.06212735   0.06426019   0.1489170
Blog        0.2259241  0.06046431   0.06217926   0.1440781

博客文章的结果等同于HC1。如果博客文章类似于你的Stata输出，Stata使用HC1。
Frank Harrel的函数产生类似于HC0的结果。据我所知，这是第一个提出的解决方案，当您查看vignette(sandwich)或?sandwich::vcovHC中提到的文章时，其他方法具有稍微更好的属性。它们在自由度调整方面有所不同。还请注意，对robcov(., method = "efron")的调用类似于HC3。

无论如何，如果您想要相同的输出，请使用HC1或适当调整方差协方差矩阵。毕竟，在查看vignette（sandwich）以了解不同版本之间的差异后，您会发现只需缩放一个常数即可从HC1转换为HC0，这应该不难。顺便说一句，需要注意的是，由于在存在影响观察值时具有更好的小样本性质和行为，因此通常首选HC3或HC4。所以，您可能想要更改Stata中的默认设置。

您可以通过将其提供给适当的函数（例如lmtest :: coeftest 或car :: linearHypothesis ）来使用这些方差协方差矩阵。例如：

library(lmtest)
coeftest(fitlm, vcov=vcovHC(fitlm, "HC1"))

t test of coefficients:

              Estimate Std. Error t value  Pr(>|t|)    
(Intercept)   1.855997   0.225924  8.2151 1.038e-13 ***
Sepal.Width   0.650837   0.060464 10.7640 < 2.2e-16 ***
Petal.Length  0.709132   0.062179 11.4046 < 2.2e-16 ***
Petal.Width  -0.556483   0.144078 -3.8624 0.0001683 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

对于集群鲁棒标准误差，您需要调整三明治的肉（请参见?sandwich），或查找一个function来执行此操作。已经有几个来源解释如何以详细到极致的细节方式进行这个操作，使用适当的代码或函数。我没有必要在这里重新发明轮子，所以我跳过这一部分。

还有一个相对较新且方便的软件包，用于计算线性模型和广义线性模型的集群稳健标准误差。请参见此处。