While I was reading @renkun-ken tutorial on pipeR, I happened to also spot this nice post from Zev Ross Recreate the GAM partial regression smooth plots from R package mgcv (with a little style).
Since I have not had a chance to use all the new features considered in pipeR 0.5, I thought recreating this post all in one long pipe would be a challenge and a great learning experience.
| # exercise in piping | |
| # using the very fine work of Zev Ross who deserves all the credit | |
| # http://zevross.com/blog/2014/09/15/recreate-the-gam-partial-regression-smooth-plots-from-r-package-mgcv-with-a-little-style/ | |
| # devtools::install_github("renkun-ken/pipeR@0.5") | |
| library(pipeR) | |
| library(dplyr) | |
| library(ggplot2) | |
| library(mgcv) | |
| "http://zevross.com/blog/wp-content/uploads/2014/08/chicago-nmmaps.csv" %>>% | |
| read.csv( as.is=T ) %>>% | |
| mutate( date = as.Date(date) ) %>>% | |
| filter( date > as.Date("1996-12-31") ) %>>% | |
| mutate( year = substring(date,1,4) ) %>>% | |
| # model - temperature against ozone | |
| # will choose to use pipeR side effect | |
| ( | |
| ~ thegam = gam(o3~s(temp), data=.) %>>% | |
| ( ~ plot(., residuals=TRUE, main="Yuck, not a nice plot") ) | |
| ) %>>% | |
| # create a sequence of temperature that spans your temperature | |
| (~ data.frame( temp = seq( min(.$temp), max(.$temp), length=300 ) ) %>>% | |
| (~t~{ | |
| # predict only the temperature term (the sum of the | |
| # term predictions and the intercept gives you the overall | |
| # prediction) | |
| predict( | |
| thegam | |
| , type="terms" | |
| , newdata=t | |
| , se.fit=TRUE | |
| ) %>>% | |
| (~ | |
| # plot the temperature smooth but leave blank for | |
| # now so that we can add the line on top of the polygon | |
| plot(t$temp, .$fit, type="n", lwd=3, xlim=c(-3,90), ylim=c(-20,30), | |
| main="Ahhh, definitely better", | |
| ylab=paste("s(temp,", round(sum(thegam$edf[-1]),2), ")", sep="")) | |
| ) %>>% | |
| (~ | |
| # For the confidence grey polygon | |
| polygon(c(t$temp, rev(t$temp)), | |
| c(.$fit+1.96*.$se.fit,rev(.$fit-1.96*.$se.fit)), col="grey", | |
| border=NA) | |
| ) %>>% | |
| (~ | |
| lines(t$temp, .$fit, lwd=2) | |
| ) | |
| }) | |
| ) %>>% | |
| (~ origdata ~( | |
| predict(thegam,type="terms") + residuals(thegam) %>>% | |
| (~ | |
| points(origdata$temp,., pch=16, col=rgb(0, 0, 1, 0.25)) | |
| ) | |
| )) %>>% | |
| # if you want to add the rug plot | |
| ( rug(.$temp) ) | |