Hi again,

As a follow-up question... If I need to calculate a risk difference with a
CI, I imagine the solution would be similar? Here's how I would get the
point estimate:

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(1,2),labels=c("cat1","cat2")),
                   outcome=c(0,0))

library(boot)
risks <- inv.logit(predict(model,newd))
risk.diff <- risks[2] - risks[1]

Many thanks,

Dominic C.

2012/3/20 Dominic Comtois <dominic.comtois@gmail.com>

> Case solved. Thanks a lot Peter!
>
> Dominic C.
>
>
> -----Message d'origine-----
> De : peter dalgaard [mailto:pdalgd@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.mes@cbs.dk  Priv: PDalgd@gmail.com
>
>

	[[alternative HTML version deleted]]

