glm和相對風險 - 在R的重複Stata代碼

Question

任何人都可以幫助我在R中複製這些相對風險計算（及其置信區間）嗎？

在Stata中使用的類似程序描述爲here。任何人都可以告訴我如何在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)) 

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 

注意的是，在第二glm()呼叫(link="log")是沒有必要的，因爲它是默認的鏈接功能family="poisson"時。

爲了推導出「穩健」的估計，我必須讀this，這非常有幫助。您必須使用夾層包中的vcovHC()函數來獲得與glm()計算的方差矩陣不同的方差 - 協方差矩陣，並使用該函數計算標準誤差和參數估計值。

「穩健」的估計值與我從Stata獲得的估計值幾乎相同，直到小數點後第三位。所有其他結果完全相同;運行代碼並親自查看。噢，還有一件事：當你使用非分層設計的glm（）找到你自己的方式時，然後映射你的方式穿越survey包，其中包含這個和其他複雜設計修改後的分析函數的對應。我還建議您閱讀Thomas Lumley的書「複雜調查：使用R分析指南」。