Generalized additive models are an established tool to model correlations for nonlinear covariates without too much hazzle with the form of the assosiation between predictor and response. It is a great straightforward tool, which leaves most of the work to the computer.
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! To achieve this one needs customization (colors, understandable axes-labels, scaling).
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:

Stratified Additive Regression Model
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")

Stratified Additive Regression Model on Response Scale
[...] default termplot() function and specifiy se=T to get the confidence limits as lines. I showed in a recent post how to plot the fitted GAM smooth and the confidence limits [...]