glm et risque relatif -répliquer le code Stata dans R

Question

Quelqu'un pourrait-il m'aider à reproduire ces calculs de risque relatif (et son intervalle de confiance) dans R? Un procédé similaire utilisé dans Stata est décrit here. Quelqu'un pourrait-il me montrer comment faire cela dans R (mes données ont des groupes et des strates mais j'ai pris un exemple plus simple)? J'ai essayé la fonction relrisk.est mais je préfère utiliser le package d'enquête car il gère des conceptions très complexes. J'aimerais aussi comparer les estimations de Stata et R. J'utilise Poisson comme suggé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

Votre exemple est un jeu de données simple, vous n'avez donc pas du tout besoin d'utiliser le package de sondage. Je suggère également que, en apprenant la régression multiple avec R, vous commencez avec des exemples plus simples et développez graduellement votre compréhension des méthodes impliquées.

Après tout, mon opinion est que les philosophies de Stata et R diffèrent. Stata jette facilement une tonne d'informations sur votre visage, sans que vous sachiez ce que cela signifie ou comment il est dérivé. D'un autre côté, R peut être tout aussi puissant (voire plus) et plus polyvalent, mais il manque l'automatisation de Stata et vous force à aller lentement pour obtenir les résultats que vous voulez.

Voici une traduction plus littérale de votre code 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 

Notez que le (link="log") dans le second glm() appel n'est pas nécessaire, car il est la fonction de lien par défaut lorsque family="poisson".

Pour dériver les estimations «robustes», j'ai dû lire this, ce qui était très utile. Vous devez utiliser la fonction vcovHC() dans le package sandwich pour obtenir une matrice de variance-covariance différente de celle calculée par glm(), et l'utiliser pour calculer les estimations d'erreur et de paramètre standard.

Les estimations «robustes» étaient presque identiques à ce que j'ai obtenu de Stata, jusqu'à la troisième décimale. Tous les autres résultats étaient complètement identiques; exécutez le code et voyez par vous-même. Oh, et encore une fois: quand vous trouvez votre chemin en utilisant glm() avec des motifs non stratifiés, faites votre chemin à travers le paquet survey, qui contient des contreparties de cette fonction et d'autres fonctions d'analyse modifiées. Je vous recommande également de lire le livre de Thomas Lumley «Complex Surveys: un guide d'analyse en utilisant R».