Hi, In fitting GAMs to assess environmental preferences, I use the part of the fit where the lower confidence interval is above zero as my criterion for positive association between the environmental variable and species abundance. However I like to plot this on the original scale of species abundance. To do so I extract the fit and SE using predict.gam.
Lately I compared more carefully the plots I obtain in this way and those obtained with plot.gam and noticed differences which I do not understand. To avoid sending a large dataset I took an example from gam Help to illustrate this. Was I wrong to believe that the fit and its confidence band should behave the same way on both scales? Thanks in advance, Denis Chabot ####################### library(mgcv) set.seed(0) n<-400 sig<-2 x0 <- runif(n, 0, 1) x1 <- runif(n, 0, 1) x2 <- runif(n, 0, 1) x3 <- runif(n, 0, 1) f0 <- function(x) 2 * sin(pi * x) f1 <- function(x) exp(2 * x) f2 <- function(x) 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10 f <- f0(x0) + f1(x1) + f2(x2) g<-exp(f/4) y<-rpois(rep(1,n),g) mean.y <- mean(y) gam2 <- gam(y~ s(x2), poisson) # to plot on the response scale val.for.pred <- data.frame(x2=seq(min(x2), max(x2), length.out=100)) pred.2.resp <- predict.gam(gam2, val.for.pred ,type="response", se.fit=TRUE) lower.band <- pred.2.resp$fit - 2*pred.2.resp$se.fit upper.band <- pred.2.resp$fit + 2*pred.2.resp$se.fit pred.2.resp <- data.frame(val.for.pred, pred.2.resp, lower.band, upper.band) # same thing on term scale pred.2.term <- predict.gam(gam2, val.for.pred ,type="terms", se.fit=TRUE) lower.band <- pred.2.term$fit - 2*pred.2.term$se.fit upper.band <- pred.2.term$fit + 2*pred.2.term$se.fit pred.2.term <- data.frame(val.for.pred, pred.2.term, lower.band, upper.band) # it is easier to compare two plots instead of looking at these two data.frames plot(gam2, residuals=T, pch=1, cex=0.7) abline(h=0) plot(y~x2, col=grey(0.5)) lines(fit~x2, col="blue", data=pred.2.resp) lines(lower.band~x2, col="red", lty=2, data=pred.2.resp) lines(upper.band~x2, col="red", lty=2, data=pred.2.resp) abline(h=mean.y,lty=5,col=grey(0.35)) ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
