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.comt...@gmail.com>
> 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
>
>
[[alternative HTML version deleted]]
______________________________________________
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.
