gl和相对风险 - 重复R中的Stata代码

Question

任何人都可以帮我复制R中的这些相对风险计算(及其置信区间)吗？

这里描述了Stata中使用的类似过程.任何人都可以告诉我如何在R中做到这一点(我的数据有集群和分层,但我采取了一个更简单的例子)？我已经尝试了relrisk.est功能,但我宁愿使用调查包,因为它处理非常复杂的设计.我还想比较Stata和R估计.我在这里建议使用泊松.

###STATA CODE
use http://www.ats.ucla.edu/stat/stata/faq/eyestudy 
tabulate carrot lenses
*same as R binomial svyglm below
xi: glm lenses carrot, fam(bin) 
*switch reference code
char carrot[omit] 1
xi: glm lenses i.carrot, fam(poisson) link(log) nolog robust eform

###R
library(foreign)
library(survey)
D<-read.dta("http://www.ats.ucla.edu/stat/stata/faq/eyestudy.dta")
table(D$lenses,D$carrot)
D$wgt<-rep(1,nrow(D))
Dd<-svydesign(id=~1,data=D,weights=~wgt)
#change category and eform....?
svyglm(lenses~carrot,Dd,family=binomial)
svyglm(lenses~carrot,Dd,family=quasipoisson(log))

Answer 1

您的示例是一个简单的数据集,因此您根本不需要使用调查包.我还建议,在使用R学习多元回归时,您可以从更简单的示例开始,逐步建立您对所涉及方法的理解.

毕竟,我的观点是Stata和R的哲学不同.Stata很容易在你面前抛出大量的信息,而你却不知道它意味着什么或它是如何衍生出来的.另一方面,R可以像(甚至更多)一样强大且更通用,但缺乏Stata的"自动化"并迫使你慢慢获得你想要的结果.

以下是Stata代码的更直译:

require(foreign)
data <- read.dta("http://www.ats.ucla.edu/stat/stata/faq/eyestudy.dta")
with(data, table(carrot, lenses))
glm.out.1 <- glm(lenses ~ carrot, family="binomial", data=data)
summary(glm.out.1)
logLik(glm.out.1)   # get the log-likelihood for the model, as well
glm.out.2 <- glm(lenses ~ factor(carrot, levels=c(1,0)), family="poisson"(link="log"), data=data)
summary(glm.out.2)
logLik(glm.out.2)
exp(cbind(coefficients(glm.out.2), confint(glm.out.2)))

# deriving robust estimates. vcovHC() is in the sandwich package.
require(sandwich)
robSE <- sqrt(diag(vcovHC(glm.out.2, type="HC")))[2]
rob <- coef(glm.out.2)[2]
rob <- exp(c(rob, rob+qnorm(0.025)*robSE, rob-qnorm(0.025)*robSE))
names(rob) <- c("", "2.5 %", "97.5 %")
rob

请注意,(link="log")第二次glm()调用中没有必要,因为它是默认的链接函数family="poisson".

为了得出"稳健"的估计,我必须阅读这个,这非常有帮助.您必须使用vcovHC()三明治包中的函数来获得与计算的方差 - 协方差矩阵不同的方差 - 协方差矩阵glm(),并使用它来计算标准误差和参数估计值.

"稳健"估计几乎与我从Stata得到的估计相同,一直到小数点后三位.所有其他结果完全相同; 运行代码并亲眼看看.

哦,还有一件事:当您使用非分层设计的glm()找到自己的方式时,然后在整个survey软件包中映射您的方式,其中包含针对复杂设计修改的此对象和其他分析函数.我还建议你阅读Thomas Lumley的书"复杂调查:使用R分析指南".