GAM optimization methods in mgcv R package - which to choose?

Question

In mgcv there are various methods to finding the smoothing parameter, lambda, such as GCV and ML/REML. GCV works by minimizing predictive error, but is subject to under/over-smoothing. ML/REML are not as reliable, but work has been done to test their overall performance in relation to GCV (Wood, 2011).

When method = "GCV.Cp" (the default) is used in gam(), gam.check automatically uses the 'magic' as its optimizer. When method = "REML", gam.check() uses bfgs optimization.

I am also looking to use Shape-Constrained Additive Models (SCAMs) which employs bfgs optimization in scam.check() and GCV in scam() as defaults. Not being an expert in GAM theory by any means, I am confused on which methods/optimizers to use in order to make comparing unconstrained GAMs to SCAMs as straightforward as possible, ensuring that models can indeed be compared successfully using AIC.

I am fitting ecological data and will be presenting my idea at a conference in November to a group of biologists (who don't know what GAMs/scam are).

Any advice on comparing GAMs/SCAMs is greatly appreciated.

Answer 1

This question is quite old now so perhaps some of the detail of {mgcv} has changed but I don't believe that is the case. That given, there are some inconsistencies in the stated facts.

{mgcv} selects the values of smoothing parameters by finding the values of the $\lambda$s and the values of $\boldsymbol{\beta}$ (by which I mean all the coefs for the parametric terms plus the coefs for the basis functions) by minimising a criterion. Which criterion is used is controlled via the method argument to GAM, and {mgcv} provides many:

• GCV.Cp (the default) which is GCV for models with unknown scale parameter and Mallows Cp/UBRE/AIC for models with known scale (poisson, binomial)
• GAV.Cp is as per above but using Generalized Approximate Cross Validation in place of GCV.
• REML or ML for REML or ML estimation,
• P-REML or P-ML, as above but using a Pearson estimate for the scale parameter

So them's the criteria we can minimise to select the smoothing parameters and other model coefficients.

But we need to have a method for actually minimising those criteria, whichever we choose. Given certain information about the model and the criterion being minimised, these optimization methods are what actually do the process of taking steps in the parameter space to find the minimum of the criterion. Which optimisation algorithm is used is controlled by the optimizer argument.

There are (were) essentially three choices, but one of these is now deprecated (the first)

• perf for performance iteration
• outer for outer iteration
• efs for the extended Fellner Schall method

The first of these was the original choice but is now deprecated as it doesn't work very well.

The outer iteration selects some value for the smoothing parameters and then via an inner iteration finds the values of the model coefficients given the selected values of the smoothing parameters. This needs quite a lot of information about the derivatives of parameters etc to allow it to update the smoothing parameters at the outer step of the next iteration of the algorithm (i.e. how to moth in the smoothing parameter space such that the criterion (method) we chose will be reduced.

This outer algorithm can use many different optimisation algorithms to find the model parameters. Which optimisation algorithm is used is specified in the second element passed to optimizer. The default is newton for Newton's method, but bfgs, optim, nlm, and nlm.fd are allowed options.

IIRC, Simon Wood has shown that under certain assumptions, the outer iteration approach is gauranteed to converge.

The final option, efs for extended Fellner Schall method is a more recent addition to {mgcv} and provides an algorithm that needs less information about derivatives of values used in the optimisation than outer. As such, it can be used with models fitted with new families without one having to tediously find and write out equations for the 3rd and 4th order derivatives of the required values/functions that outer needs. The downside is that the EFS algorithm is not guaranteed to converge.

OK, that's a lot of info and I have skipped all the detail - read Simon's GAM book for that.

I'm not sure why you think REML/ML selection is "not reliable"? I think REML/ML might have a tendency to oversmooth relative to GCV is the latter is working correctly, but GCV is known to perform badly and inconsistently in a range of situations for the problem of selecting smoothing parameters. As such, Simon and many of us in the ecology-GAM community recommend using REML or ML smoothing selection. GCV is still useful is you want to focus on minimising prediction error when choosing smoothing parameters, but beware it can undersmooth and it doesn't handle concurvity as well.

gam.check() and scam.check() are not using different optimization algorithms to gam() and scam() - all they are doing is reporting what optimization method was used (optimizer).

When gam.check() reports that it is using magic I think this is an inconsistency as {mgcv} has been iteratively developed. Basically this means GCV with outer iteration using Newton's method and hence is the name given to the combo optimizer = c("outer", "newton").

scam() has a different default; it minimises the GCV.Cp criterion (I think this only option in {scam}) using outer iteration with the bfgs optimiser: optimizer = c("outer", "bfgs"). scam() has optimizer = "bfgs" as the default, but that's just because the "outer"` bit is forced.

With method = "REML" or method = "ML" and gam(), gam.check() will actually report:

Method: REML   Optimizer: outer newton

This is the same combination of optimizer and smoothing parameter selection algorithm as the "GCV.Cp" default, but for historical reasons it is reported separately. magic is an old function in {mgcv} that was there for "outer" + "newton" with "GCV.Cp" long before REML or ML smoothing parameter selection were options in {mgcv}.

To sum up:

1. use method = "REML" unless you are comparing models (likelihood ratio tests, AIC) with different fixed effects (i.e. any model containing non-fully penalized smooths - i.e. most default smooths, exceptions are "fs" and "re" bases),
2. use method = "ML" if you do want to compare models with different fixed effects or containing non-fully penalized smooths
3. use "outer" iteration unless you have reason to change to "efs" (reasons to change might be to speed up fitting for very large data using more exotic families in {mgcv} including the location-scale-shape families)
4. use "newton" as the optimisation algorithm unless you have convergence issues. While convergence issues usually indicate that you are trying to fit too many terms for the number of data or are fitting models that contain parameters that are poorly identified, sometimes one optimisation algorithm might struggle somewhere in the parameter space while another one may not. If you get convergence issues that you can't track down to trying to fit too complex a model or with poorly identified parameters, you could try switching to one of the other algorithms and see if the issue persist - if they do it's likely that you are asking too much of the data with the model you are trying to fit.