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]]

