# How2plot nicer GAM curves

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") ```

## 2 thoughts on “How2plot nicer GAM curves”

1. Johannes says:

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")