In case someone else is interested, failed to cc to the list. Sorry
F

Dear Adam,

The invlogit only works when you have a predicted odds, invlogit on an
odds-ratio will be non-sensical, e.g. invlogit(a:b) wont do any good.
Similarly for the exp.

You are using continuous covariates, in that case observations lying nearby
will influence the estimate of what you predict which is the reason why

nn<-predict(b.glm, type="response")
mm<-predict(p.glm, type="response")/test$total
all.equal(mm,nn)

aren't totally concordant since the weighting will be different (by using
different nonlinear transformations of the probability). But they play in
the same league (and are on the correct scale):

plot(nn,mm)
abline(c(0,1), col="red")

Actually, there may be a reason that one would like to use poisson with
offset, namely in the case when the estimated probability is 1. Here the
binomial will find difficulties with SE of the estimates etc. whereas the
poisson model oughtn't. It doesn't seem to go that well compared to the
Clopper-Pearson CI, although it is an improvement in relation to the
binomial model (see code below).

I agree with the beautiful words of Spencer Graves, although I don't think
that one should flame each other. (Thank you Ben for this fortune cookie!).
Also, I don't think that one should refrain from posting ones ideas even if
they may be wrong. We are all wrong sometimes. By communicating we can
improve and get closer to the truth/good models, which is what I like so
much about the R-community.

Similarly to continuous covariates mixed models will be different if one
use a mixed-logistic or a mixed-poisson with offset. Hopefully they will be
similar enough.

Cheers,
Fredrik


invlogit<-function(x) exp(x)/(1+exp(x))
n<-10
nn<-ceiling(exp(rnorm(n)+3))
pp<-rep(c(0.25,1),each=n/2)
ca<-rbinom(rep(1,n),nn,prob=
pp)
co<-nn-ca
fa2<-gl(2,n/2)


#binomial model
test.glm<-glm(cbind(ca,co)~fa2, family=binomial)
#conditional poisson
ptest.glm<-glm(ca~fa2+offset(log(nn)), family=poisson)

summary(ptest.glm)
exp(coef(ptest.glm))
invlogit(coef(test.glm))
all.equal(exp(coef(ptest.glm)), invlogit(coef(test.glm)))
# wrong of course
# which is approximately what you do when you take invlogit(a:b) etc.
# consider instead
a<-cumsum(coef(test.glm))
b<-cumsum(coef(ptest.glm))
all.equal(invlogit(a),exp(b))

#Clopper-Pearson CI (which have a coverage probability of at least .95)
p1<-binom.test(sum(ca[1:5]), sum(nn[1:5]))$conf.int
p2<-binom.test(sum(ca[6:10]), sum(nn[6:10]))$conf.int
print(rbind(p1,p2))

#Wald CI
Vpt<-vcov(ptest.glm)
w1<-c(exp(coef(ptest.glm)[1]- sqrt(c(1,0) %*% Vpt %*% c(1,0))*2),
      exp(coef(ptest.glm)[1]+ sqrt(c(1,0) %*% Vpt %*% c(1,0))*2))
w2<-c(exp(sum(coef(ptest.glm))- sqrt(c(1,1) %*% Vpt %*% c(1,1))*2),
      exp(sum(coef(ptest.glm))+ sqrt(c(1,1) %*% Vpt %*% c(1,1))*2))
print(rbind(w1,w2))

library(multcomp)
A<-matrix(c(1,0,1,1), nrow=2,ncol=2, byrow=TRUE)
test.glht<-glht(test.glm, linfct=A)
ptest.glht<-glht(ptest.glm, linfct=A)
invlogit(confint(test.glht)$confint)
exp(confint(ptest.glht)$confint)
#not that beuatiful either

	[[alternative HTML version deleted]]


