R help - Mar 2012 - glm: getting the confidence interval for an Odds Ratio, when using predict()

Say I fit a logistic model and want to calculate an odds ratio between 2
sets of predictors. It is easy to obtain the difference in the predicted
logodds using the predict() function, and thus get a point-estimate OR. But
I can't see how to obtain the confidence interval for such an OR.

 

For example:

model <- glm(chd ~age.cat + male + lowed, family=binomial(logit))

pred1 <- predict(model, newdata=data.frame(age.cat=1,male=1,lowed=1))

pred2 <- predict(model, newdata=data.frame(age.cat=2,male=0,lowed=0))

OR <- exp(pred2-pred1) 

 

Thanks


	[[alternative HTML version deleted]]

On Mar 19, 2012, at 03:32 , Dominic Comtois wrote:
> Say I fit a logistic model and want to calculate an odds ratio between 2
> sets of predictors. It is easy to obtain the difference in the predicted
> logodds using the predict() function, and thus get a point-estimate OR. But
> I can't see how to obtain the confidence interval for such an OR.
> 
> 
> 
> For example:
> 
> model <- glm(chd ~age.cat + male + lowed, family=binomial(logit))
> 
> pred1 <- predict(model, newdata=data.frame(age.cat=1,male=1,lowed=1))
> 
> pred2 <- predict(model, newdata=data.frame(age.cat=2,male=0,lowed=0))
> 
> OR <- exp(pred2-pred1) 

There's no trivial way since you need the covariance of pred2 and pred1 to
calculate the variance of the difference.

I think you can proceed somewhat like as follows (I can't be bothered to
test it without a reproducible example to start from. You may need to throw in a
few explicit t() and as.vector() here and there.)

newd   <- data.frame(age.cat=c(1,2),male=c(1,0),lowed=c(1,0))
M <- model.matrix(model, data=newd)
V <- vcov(model)
contr <- c(-1,1) %*% M
se <- contr %*% V %*% contr 

OR.ci <- exp(pred2 - pred1 + qnorm(c(.025,.50,.975))*se)

(Sanity check:  contr %*% coef(model)  should be same as  pred2 - pred1 )

I'm not sure how general the model.matrix trick is. It works in cases like
> mm <- glm(ff, data=trees)
> model.matrix(mm, data=trees[1,])
  (Intercept) log(Height) log(Girth)
1           1    4.248495   2.116256
attr(,"assign")
[1] 0 1 2

but I see that there are cases where a "data" argument may be ignored.
If that is the case, then you may have to construct the "contr" vector
by hand.
> 
> 
> 
> Thanks
> 
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at 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.
-- 
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