[R] glm: getting the confidence interval for an Odds Ratio, when using predict()

[Dominic Comtois]

Case solved. Thanks a lot Peter!

Dominic C.


-----Message d'origine-----
De : peter dalgaard [mailto:pdalgd at gmail.com] 
Envoyé : 20 mars 2012 07:57
À : Dominic Comtois
Cc : r-help at 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.mes at cbs.dk  Priv: PDalgd at gmail.com