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…
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.