[R] glm: getting the confidence interval for an Odds Ratio, when using predict()
peter dalgaard
pdalgd at gmail.com
Tue Mar 20 12:56:31 CET 2012
[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
More information about the R-help
mailing list