Generalized additive models visualize potential non-linear associations between a predictor and a response variable.
The default plotting method produces clean plots for all covariates in the model (or a selection) but: They do not have presentation quality by any means in terms of colors, understandable axes-labels, scaling, etc.
This example shows a customization variant for a additive regression model with one covariate. The goal was to display the absolute value of the response variable on the y-axis and not the “difference from intercept” which is default.
This is only meaningful for a single covariate in the model.
1 The default plot.gam() method
MyGAM1<- with(MyData[MyData$Strata==1,], gam(Y ~ s(Covariate)))
MyGAM0<- with(MyData[MyData$Strata==0,], gam(Y ~ s(Covariate)))
par( mfcol=c(1,2))
plot(MyGAM0)
plot(MyGAM1)
This is the resulting plot:
2 The fancy way
1. Extract the values of the model response from the GAM object:
response1 <- predict(MyGAM1, type="response", se.fit=T)
response0 <- predict(MyGAM0, type="response", se.fit=T)
2. Print the response values against the covariate (note: this works just with one covariate)
par(mfcol=c(1,1))
plot(0, type="n", bty="n", main="Fancy GAM plot", xlab="MyCovariate", ylab="MyResponse", lwd=3,ylim=c(0,60), xlim=c(0,200))
legend("bottomright", bty="n", lwd=5, col=c("green","red"), legend=c("Strata = 0", "Strata = 1"))
lines(sm.spline(MyGAM1$model$Covariate , response1$fit) , lwd = 3 , col = "red")
lines(sm.spline(MyGAM1$model$Covariate , response1$fit+1.96*response1$se) , lty = 3 , lwd = 2 , col = "red")
lines(sm.spline(MyGAM1$model$Covariate , response1$fit-1.96*response1$se) , lty = 3 , lwd = 2 , col = "red")
lines(sm.spline(MyGAM0$model$Covariate , response0$fit) , lwd = 3 , col = "green")
lines(sm.spline(MyGAM0$model$Covariate, response0$fit + 1.96 * response0$se) , lty = 3 , lwd = 2, col = "green")
lines(sm.spline(MyGAM0$model$Covariate, response0$fit - 1.96 * response0$se) , lty = 3 , lwd = 2 , col = "green")
abline(h=gam.dm1$coefficients[1], lty=2, lwd=1, col="red")
abline(h=gam.dm0$coefficients[1], lty=2, lwd=1, col="green")
Rforge 
This produces really neat graphs, but it does not reproduce a plot made with plot.gam(). plotgam() plots confidence bands for the parameter estimates while this code produces confidence intervals for fitted values (which are necessarily wider). This doesn’t make it a bad solution, but one should keep it in mind when interpreting the result.
See this example code:
library(mgcv)
set.seed(0)
dat <- gamSim(1, n=400, dist="normal", scale=2)
gam1 <- gam(y~s(x0), data=dat)
plot(gam1)
newdat <- data.frame(x0=seq(from=0, to=1, by=0.01))
fits1 <- predict(gam1, newdata=newdat, type="response", se.fit=TRUE)
lines(newdat$x0, fits1$fit-mean(predict(gam1, newdata=dat)), col="red")
lines(newdat$x0, fits1$fit-mean(predict(gam1, newdata=dat)) + 1.96*fits1$se.fit, col="red")
lines(newdat$x0, fits1$fit-mean(predict(gam1, newdata=dat)) – 1.96*fits1$se.fit, col="red")