What interpretation do REML/fREML values provide in generalized additive models (GAMs)?

I'm continuing my slow trudge through Simon Wood's book on generalized additive models (GAMs), and it has given me some new useful insights. However, I am still confused after reading through Chapter 2 about what value the REML/fREML estimates provide in GAM summaries. For example, I have fit these two models using REML below:

#### Load Libraries ####
library(mgcv)
library(gamair)

data("wesdr")
wes <- tibble(wesdr)

#### Fit Candidate Models ####
fit.1 <- gam(
ret ~ s(dur, bs = "cr"),
data = wes,
method = "REML"
)

fit.2 <- gam(
ret ~ s(dur, bs = "tp"),
data = wes,
method = "REML"
)

#### Summarize ####
summary(fit.1)
summary(fit.2)


The REML values for these two GAMS are the following:

• fit.1: 467.41
• fit.2: 469.33

I rarely every hear or see anything about these values from most articles and videos I have seen on GAMs. I also see pretty much no explanation about these values elsewhere here. However, I feel like they exist for a reason. What utility do these values have?

All the fitting methods in gam()/bam() return some form of smoothness selection score:

1. GCV
2. UBRE (AIC)
3. GACV
4. REML score
5. ML score
6. ...

and this is what is displayed in the summary output.

Their main utility is as a quick way to compare models (in the same way you might compare models based on AIC); the lower the score the "better" the model. Note that for ML or REML scores, these are not the log-likelihood of the data given at the MLEs of the model parameters: you'll get slightly different values from logLik():

r$> logLik(fit.1) 'log Lik.' -457.0239 (df=6.872103) r$> logLik(fit.2)
'log Lik.' -456.9475 (df=7.106575)


and I'm not exactly sure why (probably the REML/ML scores aren't including certain other parameters that are in the log likelihood of the data.)

As to your specific example, the model using thin plate regression spline fits slightly better (it has a slightly lower negative log restricted likelihood score), but uses a bit more wiggliness than the CRS model (fit.1). But the difference is negligible and certainly not worth bothering in this case.

library("gratia")
compare_smooths(fit.1, fit.2) |>
draw()