Thanks for your answer, much appreciated.
This ain't trivial indeed. I worked my way through it, until I got a "non
conformable arguments" error when trying to calculate the new standard
error. Since I'm not following 100% what's happening, it's hard for me to
figure out what I should do next.
Here's an example with simulated data:
x1 <- rbinom(100,1,.5)
x2 <-
factor(round(runif(100,1,5)),labels=c("cat1","cat2","cat3","cat4","cat5"))
outcome <- rbinom(100,1,.2)
model <- glm(outcome~x1+x2,family=binomial(logit))
newd <-
data.frame(x1=factor(c(0,1)),x2=factor(c("cat1","cat2")),outcome=c(0,0))
M <- model.matrix(model1, data=newd)
V <- vcov(model1)
contr <- c(-1,1) %*% M
se <- contr %*% V %*% contr
OR.ci <- exp(pred2 - pred1 + qnorm(c(.025,.50,.975))*se)
Thanks for any additional hints,
Dominic C.
2012/3/19 peter dalgaard <pda...@gmail.com>
>
> 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.
>
> --
> 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.
