[R] confidence bands for a quasipoisson glm

[David Winsemius]

On Sep 11, 2010, at 3:15 PM, Maik Rehnus wrote:

> Dear all,
>
> I have a quasipoisson glm for which I need confidence bands in a  
> graphic:
>
> gm6 <- glm(num_leaves ~  b_dist_min_new, family = quasipoisson, data  
> = beva)
> summary(gm6)
>
> library('VIM')
> b_dist_min_new <- as.numeric(prepare(beva$dist_min,  
> scaling="classical", transformation="logarithm")).
>
> My first steps for the solution are following:
>
> range(b_dist_min_new)
> x <- seq(-1.496, 1.839, by=0.01)
> newdat <- data.frame(b_dist_min_new=x)
> y <- predict(gm6, newdata=newdat, type="response")
> plot(x,y, type="l", ylim=c(0,15), lty=2, xlab="Distance [scaled  
> log.]", ylab="Number of used plant", las=1)
>
> ilogit<-function(x) exp(x)/(1 + exp(x))
> logit <-function(x) log(x/(1 - x))
>
> newdat$logitpred <- predict(gm6, newdata=newdat, type="link")

I'm puzzled. You specified that model as quasipoisson and are now  
treating it as if it were  logistic? The link is going to be log(),  
nicht wahr?

> newdat$sepred <- predict(gm6, newdata=newdat, type="link",  
> se.fit=TRUE)$se.fit
> newdat$logitlower <- newdat$logitpred-1.96 * newdat$sepred
> newdat$logitupper <- newdat$logitpred+1.96 * newdat$sepred

I'm not familiar with ilogit (sounds very useful assuming it to be an  
inverse logit), but if one were taking a first stab at an inverse  
function for quasipoisson wouldn't that be exp()?

> newdat$upper <- ilogit(newdat$logitupper)
> newdat$lower <- ilogit(newdat$logitlower)
> lines(x, newdat$lower, lty=3)
> lines(x, newdat$upper, lty=3).
>
> In this way I could find a positive correlation. But my created  
> confidence bands on the graph don't touch my regression line. Could  
> it be a technical problem or is it a mistake in the calculation?
>
> I am new here and I hope you can help to solve my problem. I could  
> not find any answers for quasipoisson glm on internet.
>
> Best regards
>
> Maik
>

David Winsemius, MD
West Hartford, CT