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

Может ли кто-нибудь помочь мне воспроизвести эти расчеты относительного риска (и его доверительный интервал) в 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))

r glm stata
person Marco M    schedule 13.01.2013    source источник

Ответы (1)


arrow_upward
7
arrow_downward

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

person Theodore Lytras    schedule 13.01.2013