[R] Prediction and plottiong of several effects against response with GAM (mgcv)

Torleif Markussen Lunde torleif at stryn.net
Wed Dec 6 11:02:31 CET 2006


Hi

Using GAM (mgcv) I want plots of x against y. Using predict.gam when 
there is only one effect x works fine. However, when i want to make 
plots of several x agains y. Can this be done? I have tried using 
model.gam$model$x2 in pred.2.resp whitout making it work.

Thanks
Torleif


Below is how I have done it:


ss <- alder==1&CPUE>0
x1 <- altitude[ss]
x2 <- start[ss]
x3 <- bunndyp[ss]
x4 <- longitude[ss]
y <- CPUE[ss]
mean.y <- mean(y)

model.gam <- gam(y~ s(x2, bs="cc") + s(x1, bs="ts") + 
s(I(x3^.5),bs="ts"), poisson, gamma=1.4)

# to plot on the response scale
val.for.pred <- data.frame(x2=seq(min(x2), max(x2), length.out=222), 
x1=seq(min(x1), max(x1), length.out=222), x3=seq(min(x3), max(x3), 
length.out=222))
pred.2.resp <- predict.gam(model.gam, 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(model.gam, 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(model.gam, residuals=T, pch=1, cex=0.7)
abline(h=0)

#To plot x2 against y
plot(y~x2, col=grey(0.8), pch=19, cex=0.2)
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=coef(model.gam)[1],lty=5,col=grey(0.35))




More information about the R-help mailing list