Question

When applying gam.check in the mgcv package, R produces some residual plots and basis dimension output. Is there a way to only produce the plots and not the printed output?

library(mgcv)
set.seed(0)
dat <- gamSim(1,n=200)
b   <- gam(y~s(x0)+s(x1)+s(x2)+s(x3), data=dat)
plot(b, pages=1)
gam.check(b, pch=19, cex=.3)
Was it helpful?

Solution

There are four plots, from top left, moving down and across we have:

  1. A QQ plot of the residuals
  2. A histogram of the residuals
  3. A plot of residuals vs the linear predictor
  4. A plot of observed vs fitted values.

In the code below, I assume b contains your fitted model, as per your example. First some things we need

type <- "deviance"  ## "pearson" & "response" are other valid choices
resid <- residuals(b, type = type)
linpred <- napredict(b$na.action, b$linear.predictors)
observed.y <- napredict(b$na.action, b$y)

Note the last two lines are applying the NA handling method used when the model was fitted to the information on the linear.predictors and y, the stored copy of the response data.

The above code and that shown below is all given in the first 10 or so lines of the gam.check() source. To view this, just enter

gam.check

at the R prompt.

Each plot is produced as follows:

QQ plot

This is produced via qq.gam():

qq.gam(b, rep = 0, level = 0.9, type = type, rl.col = 2, 
       rep.col = "gray80")

Histogram of residuals

This is produced using

hist(resid, xlab = "Residuals", main = "Histogram of residuals")

Residuals vs linear predictor

This is produced using

plot(linpred, resid, main = "Resids vs. linear pred.", 
     xlab = "linear predictor", ylab = "residuals")

Observed vs fitted values

This is produced using

plot(fitted(b), observed.y, xlab = "Fitted Values", 
     ylab = "Response", main = "Response vs. Fitted Values")

OTHER TIPS

There are now the two packages gratia and mgcViz which have functions to produce the gam.check output as ggplots which you can store as an object. The former doesn't print anything to console, the latter does.

require(gratia)

appraise(b)
require(mgcViz)

b = getViz(b)
check(b)
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top