I apologize in advance for my very basic question. I would like to plot prediction intervals for my GAM model, however, instead of using a Gamma GAM, I assumed a Gaussian distribution (I standardized my response variable, which resulted in negative values). To create the prediction interval for my model, I simply changed the code below: yr <- matrix(rgamma(fv*0,shape=1/b$scale,scale=fv*scale),nrow(fv),ncol(fv)) To: yr <- matrix(rnorm(fv*0,mean=mean(data$Y),sd=sd(data$Y)),nrow(fv),ncol(fv)) # where Y is my response variable However, the prediction interval created does not capture the trend in my data (see attached figure). I imagine I also should have changed: fv <- exp(lp) where Xp <- predict(b,newdata=data.frame(x=xp),type="lpmatrix") lp <- Xp%*%br Any suggestions would be greatly appreciated! Best regards, Zofia Stagiaire Postdoctoral Universit? de Montr?al D?partement de sciences biologiques Pavillon Marie-Victorin 90, av. Vincent-d'Indy Outremont, QC H2V 2S9