Combined deviance and quantile residuals check plots for GAMs

better_check.R

library(gridBase)
library(grid)
library(mgcv)
library(statmod)
## nicer version of gam.check/rqgam.chack
# - deviance residuals for the Q-Q and histogram
# - RQR for resids vs linear pred
# - response vs fitted
# hist.p gives the quantiles of the residuals to show
better_check <- function(b, k.sample = 5000, k.rep = 200, rep = 0, level = 0.9,
    rl.col = 2, rep.col = "gray80", hist.p=c(0.01, 0.099), ...){

  # layout stuff
  opar <- par(mfrow=c(2,2))

  # grab the randomised quantile residuals
  # requires statmod package

  # need to do the right thing for mgcv's Tweedie
  if(grepl("^Tweedie",b$family$family)){
    if(is.null(environment(b$family$variance)$p)){
      p.val <- b$family$getTheta(TRUE)
      environment(b$family$variance)$p <- p.val
    }
    qres <- qres.tweedie(b)
  # and for negbin
  }else if(grepl("^Negative Binomial",b$family$family)){
    # need to set $theta
    if("extended.family" %in% class(b$family)){
      # for SNW's extended family, need to set TRUE in getTheta as theta
      # is on the wrong scale
      b$theta <- b$family$getTheta(TRUE)
    }else{
      b$theta <- b$family$getTheta()
    }
    qres <- qres.nbinom(b)
  }else{
    stop("Only negative binomial and Tweedie response distributions are supported.")
  }

  # values of the linear predictor
  linpred <- napredict(b$na.action, b$linear.predictors)

  ## get the deviance residuals
  resid <- residuals(b, type = "deviance")

  ## Q-Q plot of deviance resids
  qq.gam(b, rep = rep, level = level, type = "deviance", rl.col = rl.col,
        rep.col = rep.col, main="Quantile-quatile plot", ...)

  ## Q resids vs. linear pred
  plot(linpred, qres, main="Resids vs. linear pred.",
       xlab="linear predictor", ylab="Randomised quantile residuals",
       pch=19, cex=0.5, col=rgb(0, 0, 0, 0.6), ...)
  lines(lowess(linpred, qres), col = 2)

  ## histogram of deviance resids

  vps <- baseViewports()
  pushViewport(vps$figure)

  # partial histogram
  hist(resid[resid >= quantile(resid, hist.p[1]) &
             resid <= quantile(resid, hist.p)[2]],
       xlab = "deviance residuals", main = "Histogram of residuals",...)
  box()

  # setup the smaller plot
  pushViewport(viewport(x=-.4, y=.45 ,width=.3, height=.3))
  oopar <- par(plt = gridFIG(), new=TRUE)

  # histogram of full data
  hist(resid, xlab = "", main="", axes=FALSE, ylab="", ...)
  axis(1, cex.axis=0.5, mgp=c(3, 0.1, 0), tck=-0.02)

  # reset grid options for next plot
  popViewport(2)
  par(oopar)

  ## Response vs. Fitted Values
  plot(fitted(b), b$y,
       main="Response vs. Fitted Values",
       xlab="Fitted Values", ylab="Response",
       pch=19, cex=0.5, col=rgb(0, 0, 0, 0.6), ...)
  lines(lowess(b$fitted.values, b$y), col = 2)

  par(opar)

}

Nice! A couple of comments:

If you do opar <- par() and then at the end do par(opar) as here and the function bails out for whatever reason before the par(opar) line, you'll leave the user's graphics device in a state it wasn't in before the function call. Better to do:

op <- par(mfrow = c(2,2))
on.exit(par(op))

and don't have a call to par as the last line of the function.

To be honest it's probably easier to do

layout(matrix(1:4))
on.exit(layout(1))

But that may just be me...

For some weird reason, the Response vs Fitted (lower right) is smaller than all the other figures?