glm и относительный риск-репликация кода Stata в R

Может ли кто-нибудь помочь мне воспроизвести эти расчеты относительного риска (и его доверительный интервал) в R?

Аналогичная процедура, используемая в Stata, описаназдесь . Может ли кто-нибудь показать мне, как это сделать в R (мои данные имеют кластеры и слои, но я взял более простой пример)? Я уже пробовал релриск.функция 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))
1 5
glm

1 ответ:

Ваш пример-это простой набор данных, поэтому вам не нужно использовать пакет опроса вообще. Я бы также предложил, чтобы при изучении множественной регрессии с помощью R вы начинали с более простых примеров и постепенно наращивали свое понимание используемых методов.

В конце концов, мое мнение состоит в том, что философии Stata и R различаются. Стата легко бросает тонну информации в ваше лицо, без того, чтобы вы знали, что это означает или как оно получено. С другой стороны 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, который содержит аналоги этого и другие функции анализа, модифицированные для сложных конструкций. Я также рекомендую вам прочитать книгу Томаса Ламли "комплексные опросы: руководство по анализу с использованием R".

6