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.
Thanks. I had great fun with this, getting a lot of impressive plots. I think there is a typo in step 7, line 01; there should be a closing “)”
Absolutely. Should be correct now.
Thanks for posting it.
I really appreciate the fact that you took the time to make these kinds of plots and walk people through them…I have been trying to figure it out for a few days now, and finally, this works for my data.
I did wonder if you could help me, though…I did get the line graphed against the response variable in the y-axis, but the smoothing effect is gone, and I just have a pointy line that is graphed. Any ideas?
Thanks!
Mmh, that is rather odd, since plotHR should not plot any pointy line, but solid lines and a polygon shade for the 95% confidence limits…
You might try
termplot( coxph( Your_Formula ) , se = TRUE , rug = TRUE )instead of
plotHR( coxph( Your_Formula ) )termplot() is the default R plotting method for the functional form of the covariates… I just gave it another look and the maintainers have fixed the bug which prevented plotting of a univariate model. That was the main reason I started out with plotHR()… and it renders plotHR more or less needless …. although I do prefer the default colors and CI shades of plotHR()… and the x-axis rug customization with 1Q, median, 3Q…
You might try
termplot( coxph( YourFormula ) , se = TRUE , rug = TRUE , axes = FALSE , bty = "n" )axis( side = 1 )
axis(side = 2 , at = axTicks(2) , label = round( exp(axTicks(2)) , digits = 1) )
abline( h = 0 , lty = 3 )
If that works, your model is allright and there might be a problem with plotHR()… in which case I would gladly have some feedback.
Good luck,
Reinhard
Excellent instructions! You’re sample code, explanations, and example graph really made this easy to understand. Thank you very much!
Very helpful code! Thanks!
That’s so great! It helped me a lot! Thank you so much.