有人能帮我在R中复制这些相对风险计算(及其置信区间)吗?
here中描述了Stata中使用的类似过程。谁能告诉我如何在R中做到这一点(我的数据有簇和层,但我举了一个更简单的例子)?我尝试过relrisk.est函数,但我更喜欢使用调查包,因为它可以处理非常复杂的设计。我还想比较Stata和R估计,我使用泊松作为建议的here。
###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))发布于 2013-01-13 19:43:25
您的示例是一个简单的数据集,因此您根本不需要使用调查包。我还建议,在使用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请注意,第二个glm()调用中的(link="log")不是必需的,因为它是family="poisson"时的默认链接函数。
为了得出“稳健”的估计,我不得不阅读this,这非常有帮助。您必须使用夹心包中的vcovHC()函数来获得与glm()计算的方差-协方差矩阵不同的方差-协方差矩阵,并使用该矩阵计算标准误差和参数估计。
“稳健”的估计与我从Stata得到的几乎相同,一直到小数点后第三位。所有其他结果都是完全相同的;请运行代码并亲自查看。
哦,还有一件事:当您在非分层设计中使用glm()时,然后遍历survey包,该包包含针对复杂设计而修改的此分析函数和其他分析函数的副本。我还推荐你阅读Thomas Lumley的书“复杂调查:使用R进行分析的指南”。
https://stackoverflow.com/questions/14301999
复制相似问题