The lightblue shade denoting the 95% pointwise confidence limits of the GAM estimate is a *polygon()* object in R.

Usually one would plot the GAM model with the 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 manually.

The same approach can be used to have the confidence limits as a shade. In order to achieve this:

- Fit the model
- get the estimated values of the fitted smoothing spline
- and pointwise standard errors
- calculate the values of the upper and lower 95%-confidence limits
- plot an empty coordinate system with right scaling, labels , titles and so on. I prefer to include axes = FALSE and add the axes manually. See here for an example.
- No it gets a bit tricky with sorting the coordinates in the right order. R provides a very effective way to achieve this, though it is not easily understandable at once. We create two indices on the dataset
**i.for**and**i.back**which number the dataset according to the independend variable on the x-axis – forward and backward respectively: - The points of the shade polygon follow here first the upper confidence line (ucl) forward and then the lower confidence line backward. The coordinates of the shade-polygon are
- First plot the polygon with the confidence shade – otherwise the GAM plot vanishes behind it.
- now the mainplot
**afterwards** - add a horizontal line marking the intercept ( the mean of the covariate) – note: univariate only. I will post a remark on this later…

`library(mgcv)`

model <- gam( MyResponse ~ MyCovariate , data = MyDataset)

`fit <- predict( model , se = TRUE )$fit`

`se <- predict( model , se = TRUE)$se.fit`

`lcl <- fit - 1.96 * se`

ucl <- fit + 1.96 * se

`plot( 0 , type = "n" , bty = "n" , xlab = "The x-axis lable" , ylab = "The y-axis lable" , main = "The Title" , xlim = c(0,5) , ylim = c(0,3) )`

`i.for <- order( MyDataset$MyCovariate )`

i.back <- order( MyDataset$MyCovariate , decreasing = TRUE )

`x.polygon <- c( MyDataset$MyCovariate[i.for] , MyDataset$MyCovariate[i.back] )`

y.polygon <- c( ucl[i.for] , lcl[i.back] )

`polygon( x.polygon , y.polygon , col = "#A6CEE3" , border = NA )`

`lines( MyDataset$MyCovariate[i.for] , fit[i.for], col = "#1F78B4" , lwd = 3 )`

`abline( h = mean(MyDataset$MyCovariate) , lty = 2 )`

Hope no typos. Comment if you find some.