[R] confidence bands for a quasipoisson glm

[Gavin Simpson]

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
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%