1 C02 modelling

  1. Load packages
library(qgam); library(mgcViz); library(gamair); 
data(co2s)
  1. Plot data
with(co2s, plot(c.month, co2, type="l"))

  1. Fit model for median
b <- qgam(co2~s(c.month, bs="cr", k=100), data=co2s, qu = 0.5, err = 0.1)
## Estimating learning rate. Each dot corresponds to a loss evaluation. 
## qu = 0.5..........done
  1. Use it for prediction
co2plot <- function(co2s,b) {
  fv <- predict(b,data.frame(c.month=1:543,month=c(rep(1:12,45),1:3)),se=TRUE)
  ul <- fv$fit + 2*fv$se
  ll <- fv$fit - 2*fv$se
  with(co2s,plot(c.month,co2,pch=19,cex=.3,col=2,
                 ylim=range(c(ul,ll)),xlim=c(0,550)))
  lines(1:543,fv$fit)
  lines(1:543,ul,lty=2)
  lines(1:543,ll,lty=2)
}

co2plot(co2s,b) ## nonsense predictions - extrapolation artefact

  1. Fit a better model
b1 <- qgam(co2~s(c.month,bs="cr",k=50)+s(month,bs="cc"),data=co2s, qu = 0.5,
           argGam = list(knots=list(month=c(1,13))), err = 0.1)
## Estimating learning rate. Each dot corresponds to a loss evaluation. 
## qu = 0.5.........done
  1. Predict again
co2plot(co2s,b1)

This is much better, as the short-term seasonal effect has been separated from long terms smooth terms, allowing longer range extrapolation of slow long range trend.

  1. mgcViz plotting and model checking

We convert to an gamViz object. We need to set nsim = 0 because we don’t know how to simulate from a gam model:

b1 <- getViz(b1, nsim = 0)
print(plot(b1), pages = 1) # Plotting smooth effects

qq.gam(b1)

Residuals look fairly normal, but this is only because we are fitting quantile 0.5. In fact, if we fit quantile 0.8:

b2 <- qgam(co2~s(c.month,bs="cr",k=50)+s(month,bs="cc"),data=co2s, qu = 0.8,
           argGam = list(knots=list(month=c(1,13))), err = 0.1)
## Estimating learning rate. Each dot corresponds to a loss evaluation. 
## qu = 0.8.............done
qq.gam( getViz(b2, nsim = 0) )