Case solved. Thanks a lot Peter! Dominic C.
-----Message d'origine----- De : peter dalgaard [mailto:pda...@gmail.com] Envoyé : 20 mars 2012 07:57 À : Dominic Comtois Cc : r-help@r-project.org help Objet : Re: [R] glm: getting the confidence interval for an Odds Ratio, when using predict() [Oops, forgot cc. to list] On Mar 20, 2012, at 04:40 , Dominic Comtois wrote: > I apologize for the errors in the previous code. Here is a reworked example. It works, but I suspect problems in the se calculation. I changed, from the 1st prediction to the 2nd only one covariate, so that the OR's CI should be equal to the exponentiated variable's coefficient and ci. And we get something different: Yep. Classical rookie mistake: Forgot to take sqrt() in the se. I then get > se <- sqrt(contr %*% V %*% t(contr)) > > # display the CI > exp(contr %*% coef(model) + qnorm(c(.025,.50,.975))*se) [1] 0.655531 1.686115 4.336918 > > # the point estimate is ok, as verified with > exp(model$coefficients[3]) x2cat2 1.686115 > > # however I we'd expect to find upper and lower bound equal # to the > exponentiated x2cat coefficient CI exp(confint(model))[3,] Waiting for profiling to be done... 2.5 % 97.5 % 0.6589485 4.4331058 which is as close as you can expect since the confint method is a bit more advanced than +/-2SE. -pd > x1 <- factor(rbinom(100,1,.5),levels=c(0,1)) > x2 <- > factor(round(runif(100,1,2)),levels=c(1,2),labels=c("cat1","cat2")) > outcome <- rbinom(100,1,.2) > > model <- glm(outcome~x1+x2,family=binomial(logit)) > newd <- data.frame(x1=factor(c(0,0),levels=c(0,1)), > x2=factor(c("cat1","cat2"),levels=c("cat1","cat2")), > outcome=c(1,1)) > > M <- model.matrix(formula(model), data=newd) V <- vcov(model) contr <- > c(-1,1) %*% M se <- sqrt(contr %*% V %*% t(contr)) > > # display the CI > exp(contr %*% coef(model) + qnorm(c(.025,.50,.975))*se) > > # the point estimate is ok, as verified with > exp(model$coefficients[3]) > > # however I we'd expect to find upper and lower bound equal # to the > exponentiated x2cat coefficient CI exp(confint(model))[3,] > > Many thanks, > > Dominic C. -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd....@cbs.dk Priv: pda...@gmail.com ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
