Proportion deviance explained by a predictor in a GAM fitted using mgcv with 'select=TRUE' (variable selection using double penalty)

Question

Objective: To calculate proportion deviance explained by a predictor.

Approach: Following this post by Simon Wood, the deviance explained by a predictor x1 is the difference in the deviances explained by models that (A) include and (B) do not include x1. Prof. Wood further explains it is necessary to ensure that the reduced model uses the same smoothing parameters as the full model.

Reproducible example:

library(tidyverse)
library(mgcv)

get_dev_explained = function(mod, pred) {
  # Drop 'pred' from formula
  form_wout_pred = Reduce(paste, deparse(formula(mod))) |>
    str_replace(paste0(", ", pred, "|", pred, ", "), "")
  # Drop smoothing parameter corresponding to 'pred'
  sp_wout_pred = mod$sp[str_detect(labels(terms(mod)), pred, negate = T)]
  # Fit model without predictor 'pred'
  cat("Fitting model:", form_wout_pred, "\n")
  cat("using smoothing parameters: ", paste(names(sp_wout_pred), collapse = " "), "\n")
  mod_wout_pred = update(mod, formula. = form_wout_pred, sp = sp_wout_pred)
  # Calculate % deviance explained by 'pred'
  tibble(
    pred = pred,
    with_pred = summary(mod)$dev.expl * 100,
    wout_pred = summary(mod_wout_pred)$dev.expl * 100,
    by_pred = with_pred - wout_pred)
}

set.seed(0)
n = 400
data = tibble(
  x1 = runif(n),
  x2 = runif(n),
  x3 = runif(n),
  y = x1 + x2 + x3 + rnorm(n))

mod = gam(y ~ te(x1, x2, x3), data = data, select = F)
print(mod$sp) # print smoothing parameters

te(x1,x2,x3)1 te(x1,x2,x3)2 te(x1,x2,x3)3 
1.403968e+11  4.287090e-01  8.989312e+00

There is one smoothing parameter per marginal smooth as expected.

bind_rows(
  get_dev_explained(mod, "x1"),
  get_dev_explained(mod, "x2"),
  get_dev_explained(mod, "x3")) 

Fitting model: y ~ te(x2, x3) 
using smoothing parameters:  te(x1,x2,x3)2 te(x1,x2,x3)3 
Fitting model: y ~ te(x1, x3) 
using smoothing parameters:  te(x1,x2,x3)1 te(x1,x2,x3)3 
Fitting model: y ~ te(x1, x2) 
using smoothing parameters:  te(x1,x2,x3)1 te(x1,x2,x3)2

% deviance explained by each predictor:

  pred  with_pred wout_pred by_pred
1 x1         28.6      19.9    8.67
2 x2         28.6      13.6   15.0 
3 x3         28.6      17.0   11.7

The problem: If we set select = TRUE in order to perform variable selection (the 'double penalty approach' of Marra & Wood (2011)), an extra null-space penalty term is added, effectively allowing variables to be shrunk out of the model.

mod_sel = gam(y ~ te(x1, x2, x3), data = data, select = T)
print(mod_sel$sp) # print smoothing parameters

te(x1,x2,x3)1 te(x1,x2,x3)2 te(x1,x2,x3)3 te(x1,x2,x3)4 
9.688189e+09  4.229465e-01  9.057475e+00  7.775367e-02

Question: What do the 4 smoothing parameters correspond to? i.e., which of these correspond to the marginal smooths? [Knowing this is necessary in order to use the same smoothing parameters when fitting the reduced models.]

References:

Marra, G., and S. N. Wood. 2011. Practical variable selection for generalized additive models. Computational Statistics & Data Analysis 55:2372–2387.

Answer 1

The final penalty penalizes the parametric part of the model, but in this case that includes terms x1,x2,x3,x1x2,x1x3,x2x3. The equivalent penalty will have a different meaning in your reduced models (e.g. penalizing x1,x2,x1x2).

Actually that's a generic problem with the above approach. The penalty on x1 is different between te(x1,x2) and te(x1,x2,x3) - one is over a 2D surface, the other over a 3D surface - hence the smoothing parameters are not comparable, and `holding them constant' between model fits is not meaningful.

I think that the best you can do here is to re-estimate smoothing parameter for each model fit.

Simon (mgcv author)