R help - Sep 2010 - confidence bands for a quasipoisson glm

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")   
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
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



  
-- 
GMX DSL SOMMER-SPECIAL: Surf & Phone Flat 16.000 f?r nur 19,99 Euro/mtl.!*

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

Is this something you want to have  (based on a simulated dataset)?

counts <- c(18,17,15,20,10,20,25,13,12)
#risk <- round(rexp(9,0.5),3)
risk<- c(2.242, 0.113, 1.480, 0.913, 5.795, 0.170, 0.846, 5.240, 0.648)
gm <- glm(counts ~ risk, family=quasipoisson)
summary(gm)
new.risk=seq(min(risk), max(risk),0.1)
new.risk
y <- predict.glm(gm, newdata=data.frame(risk=new.risk), se.fit=TRUE,
type="response") 
upper=y$fit+1.96*y$se.fit
lower=y$fit-1.96*y$se.fit
plot(new.risk,y$fit, type="l", col=4, lty=1, lwd=2)
lines(new.risk, upper, type="l", col=2, lty=2, lwd=2)
lines(new.risk, lower, type="l", col=2, lty=2, lwd=2)



-- 
View this message in context:
http://r.789695.n4.nabble.com/confidence-bands-for-a-quasipoisson-glm-tp2535883p2535952.html
Sent from the R help mailing list archive at Nabble.com.

On Sat, 2010-09-11 at 14:41 -0700, Peng, C wrote:
> Is this something you want to have  (based on a simulated dataset)?
> 
> counts <- c(18,17,15,20,10,20,25,13,12)
> #risk <- round(rexp(9,0.5),3)
> risk<- c(2.242, 0.113, 1.480, 0.913, 5.795, 0.170, 0.846, 5.240, 0.648)
> gm <- glm(counts ~ risk, family=quasipoisson)
> summary(gm)
> new.risk=seq(min(risk), max(risk),0.1)
> new.risk
> y <- predict.glm(gm, newdata=data.frame(risk=new.risk), se.fit=TRUE,
> type="response") 
I think you should be doing this bit on the scale of the link function,
not the response, and then transform.

y <- predict.glm(gm, newdata=data.frame(risk=new.risk), se.fit=TRUE,
                 type="link")
> upper=y$fit+1.96*y$se.fit
> lower=y$fit-1.96*y$se.fit
upper <- with(y, exp(fit + (2 * se.fit)))
lower <- with(y, exp(fit - (2 * se.fit)))
fit <- with(y, exp(fit))

plot(new.risk, fit, type="l", col=4, lty=1, lwd=2,
     ylim = range(c(upper, lower)))
lines(new.risk, upper, type="l", col=2, lty=2, lwd=2)
lines(new.risk, lower, type="l", col=2, lty=2, lwd=2)
> plot(new.risk,y$fit, type="l", col=4, lty=1, lwd=2)
> lines(new.risk, upper, type="l", col=2, lty=2, lwd=2)
> lines(new.risk, lower, type="l", col=2, lty=2, lwd=2)
HTH

G

-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
 Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%