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:
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 <- 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.
[[alternative HTML version deleted]]
______________________________________________
[email protected] 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.
