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
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
