The qq.gamViz method

Consider the following model with binomial responses:

set.seed(0)
n.samp <- 400
dat <- gamSim(1,n = n.samp, dist = "binary", scale = .33)
p <- binomial()$linkinv(dat$f) ## binomial p
n <- sample(c(1, 3), n.samp, replace = TRUE) ## binomial n
dat$y <- rbinom(n, n, p)
dat$n <- n
lr.fit <- gam(y/n ~ s(x0) + s(x1) + s(x2) + s(x3), 
              family = binomial, data = dat,
              weights = n, method = "REML")
lr.fit <- getViz(lr.fit)

We can get a QQ-plot of the residuals as follows:

qq(lr.fit, method = "simul1", a.qqpoi = list("shape" = 1), a.ablin = list("linetype" = 2))

Here method determines the method used to compute the QQ-plot, while the arguments starting with a. are lists that will be passed directly to the corresponding ggplot2 layer (geom_point and geom_abline here). We can remove the confidence intervals and show all simulated (model-based) QQ-curves as follows:

qq(lr.fit, rep = 20, showReps = T, CI = "none", a.qqpoi = list("shape" = 19), a.replin = list("alpha" = 0.2))

Importantly, mgcViz::qq.gam can handle large datasets by discretizing the QQ-plot before plotting. For instance, let’s increase n.samp in the previous example:

set.seed(0)
n.samp <- 20000
dat <- gamSim(1,n=n.samp,dist="binary",scale=.33)
p <- binomial()$linkinv(dat$f) ## binomial p
n <- sample(c(1,3),n.samp,replace=TRUE) ## binomial n
dat$y <- rbinom(n,n,p)
dat$n <- n
lr.fit <- bam(y/n ~ s(x0) + s(x1) + s(x2) + s(x3)
              , family = binomial, data = dat,
              weights = n, method = "fREML", discrete = TRUE)
lr.fit <- getViz(lr.fit)

Here the discrete argument determines whether the QQ-plot is discretized or not. Notice that we can compute the QQ-plot, store it in o and then plot it (via print.qqGam).

o <- qq(lr.fit, rep = 10, method = "simul1", CI = "normal", showReps = TRUE, 
        a.replin = list(alpha = 0.1), discrete = TRUE)
o 

The coarseness of the discretization grid is determined by the ngr argument:

o <- qq(lr.fit, rep = 10, method = "simul1", CI = "normal", showReps = TRUE,
        ngr = 1e2, a.replin = list(alpha = 0.1), a.qqpoi = list(shape = 19))
o 

We can zoom into particular areas of the plot using the zoom generic. Also, given that qq.gamViz and zoom.qqGam output objects of class plotSmooth, we can arrange them using the gridPrint:

gridPrint(o, zoom(o, xlim = c(2, 2.5), ylim = c(2, 2.5)), ncol = 2)

The QQ-plot can also be manipulated interactively using the shine generic, which transforms it into a shiny app:

# Cannot run this when building the pdf for this vignette, but do try it!
shine(o)

The check.gamViz method

The check.gamViz method is similar to mgcv::gam.check, with the difference that it produces a sequence of ggplot objects and that it sub-samples the residuals to avoid over-plotting (or stalling entirely) when dealing with large data sets. Here is an example:

set.seed(0)
dat <- gamSim(1, n = 200)

## Gu & Wahba 4 term additive model

b <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat)
b <- getViz(b)

check(b,
      a.qq = list(method = "tnorm", 
                  a.cipoly = list(fill = "light blue")), 
      a.respoi = list(size = 0.5), 
      a.hist = list(bins = 10))

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 8 iterations.
## The RMS GCV score gradient at convergence was 1.072609e-05 .
## The Hessian was positive definite.
## Model rank =  37 / 37 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##         k'  edf k-index p-value
## s(x0) 9.00 2.32    1.00    0.53
## s(x1) 9.00 2.31    0.97    0.30
## s(x2) 9.00 7.65    0.96    0.27
## s(x3) 9.00 1.23    1.04    0.67

The a.qq argument is a list that gets passed directly to qq.gamViz. Similarly, a.repoi is passed to ggplot2::geom_points and a.hist to ggplot2::geom_hist.

One dimensional checks using check1D

This function allows to verify how the residuals vary along one covariate. Consider the following model:

set.seed(4124)
n <- 5e3
x <- rnorm(n); y <- rnorm(n);
z <- as.factor( sample(letters[1:6], n, replace = TRUE) )

ob <- (x)^2 + (y)^2 + (0.2*abs(x) + 1)  * rnorm(n)
b <- gam(ob ~ s(x) + s(y) + z)
b <- getViz(b)

Here the responses variance varies a lot along x. Assume that we didn’t know this, but that we wanted to find out whether the residuals are heteroscedastic. We can start by doing the following:

ck1 <- check1D(b, "x")
ck2 <- check1D(b, "z")
gridPrint(ck1, ck2, ncol = 2)

This produces two views along x and z, but as you can see the plots are initially empty. We might want to add a layer showing the conditional distribution of the residuals along x and z and another containing a rug:

gridPrint(ck1 + l_dens(type = "cond", alpha = 0.8) + l_rug(alpha = 0.2), 
          ck2 + l_points() + l_rug(alpha = 0.2), layout_matrix = matrix(c(1, 1, 1, 2, 2), 1, 5))

The left plot suggests that the variance of the residuals might be lower in the middle (x = 0), but it is not entirely clear. The l_densCheck layer gives a more clear answer in this case:

ck1 + l_densCheck()

## It is safer to use l_densCheck when residual type is set to "tnormal" or "tunif" in check1D. 
##  Otherwise it is possible to supply a customized residuals function dFun(). See ?l_densCheck

This layers adds an heatmap proportional to {p(r|x)^1/2 − p_m(r|x)^1/2}^1/3, where r are the residuals, while p and p_m are their empirical and theoretical (model based) densities. In particular, p is estimated using the the fast k.d.e. method of Wand (1994) (implemented by the kernSmooth package) and p_m is a standard normal density here. This plot makes clear that the residuals are over-dispersed when x is far from zero.

The l_gridCheck1D provides another way of finding residuals patterns. For instance:

b <- getViz(b, nsim = 50)
gridPrint(check1D(b, "x") + l_gridCheck1D(gridFun = sd, showReps = TRUE), 
          check1D(b, "z") + l_gridCheck1D(gridFun = sd, showReps = TRUE), ncol = 2)

Here we converted b again using getViz with nsim = 50. This is because l_gridCheck1D needs some simulations to compute the confidence intervals. The simulations are done by getViz and then stored inside b. l_gridCheck1D simply bins the residuals according to their x values, and evaluates a user-defined function (sd here) over the observed and simulated residuals.

Two dimensional checks using check2D

check2D is quite similar to check1D, but looks at the residuals along two covariates. Here is an example where the mean effect follows the Rosenbrock function:

set.seed(566)
n <- 5e3
X <- data.frame("x1"=rnorm(n, 0.5, 0.5), "x2"=rnorm(n, 1.5, 1), 
                "fac"=as.factor( sample(letters[1:6], n, replace = TRUE) ))
X$y <- (1-X$x1)^2 + 100*(X$x2 - X$x1^2)^2 + rnorm(n, 0, 2)
b <- gam(y ~ te(x1, x2, k = 5), data = X)
b <- getViz(b, nsim = 50)

We start by generating two 2D views:

ck1 <- check2D(b, x1 = "x1", x2 = "x2")
ck2 <- check2D(b, x1 = X$fac, x2 = "x2") + labs(x = "fac")

Then we add the l_gridCheck2D layer:

ck1 + l_gridCheck2D(gridFun = mean)

## Warning: Computation failed in `stat_summary_hex()`.
## Caused by error in `compute_group()`:
## ! The package "hexbin" is required for `stat_summary_hex()`.

ck2 + l_gridCheck2D(gridFun = mean)

l_gridCheck2D bins the observed and simulated residuals, summarizes them using a scalar-valued function (mean here), and adds a heatmap of the observed summary in each cell, normalized using the nsim summaries obtained using the simulations. In the first plot above, the pattern in the residual means is not very well visible, due to outliers on the far right. The pattern is made more visible by zooming on the center of the distribution and by changing the size of the bins:

ck1 + l_gridCheck2D(bw = c(0.05, 0.1)) + xlim(-1, 1) + ylim(0, 3)

## Warning: Computation failed in `stat_summary_hex()`.
## Caused by error in `compute_group()`:
## ! The package "hexbin" is required for `stat_summary_hex()`.

As for smooth effect plots, we can list the available layers by doing:

listLayers( ck1 ) 

## [1] "l_dens2D"       "l_glyphs2D"     "l_gridCheck2D"  "l_points"      
## [5] "l_rug"          "l_gridQCheck2D"

The most sophisticated layer is probably l_glyphs2D which we illustrate here using an heteroscedastic model:

set.seed(4124)
n <- 5e3
dat <- data.frame("x1" = rnorm(n), "x2" = rnorm(n))
dat$y <- (dat$x1)^2 + (dat$x2)^2 + (1*abs(dat$x1) + 1)  * rnorm(n)
b <- gam(y ~ s(x1) + s(x2), data = dat)
b <- getViz(b)

ck <- check2D(b, x1 = "x1", x2 = "x2", type = "tnormal")

Similarly to l_gridCheck2D, l_glyphs2D bins the residuals according to two covariates, but the user-defined function used to summarize the residuals in each bin has to return a data.frame, rather than a scalar. Here is an example:

glyFun <- function(.d){
  .r <- .d$z
  .qq <- as.data.frame( density(.r)[c("x", "y")], n = 100 )
  .qq$colour <- rep(ifelse(length(.r)>50, "black", "red"), nrow(.qq))
  return( .qq )
}

ck + l_glyphs2D(glyFun = glyFun, ggLay = "geom_path", n = c(8, 8),
                mapping = aes(x=gx, y=gy, group = gid, colour = I(colour)), 
                height=1.5, width = 1)

Each glyph represend a kernel density of the residuals, with colours indicating whether we have more (black) or less (red) that 50 observations in that bin. It is clear that the residuals are much less variable for x ≈ 0 than elsewhere. We can do the same using binned worm-plots:

glyFun <- function(.d){
  n <- nrow(.d)
  px <- qnorm( (1:n - 0.5)/(n) )
  py <- sort( .d$z )
  clr <- if(n > 50) { "black" } else { "red" }
  clr <- rep(clr, n)
  return( data.frame("x" = px, "y" = py - px, "colour" = clr))
}

ck + l_glyphs2D(glyFun = glyFun, ggLay = "geom_point", n = c(10, 10),
                mapping = aes(x=gx, y=gy, group = gid, colour = I(colour)),
                height=2, width = 1, size = 0.2) 

Notice that worm-plots (Buuren and Fredriks, 2001) are simply rotated QQ-plots. An horizontal plot indicates well specified residual model. An increasing (decreasing) worm indicates over (under) dispersion.