Bayesian views of generalized additive modelling

TL;DR

This paper explores the connections and differences between Bayesian and frequentist approaches to generalized additive models (GAMs), emphasizing their practical applications in ecological modeling and highlighting Bayesian methods within the exttt{mgcv} framework.

Contribution

It clarifies the links between Bayesian and frequentist GAMs, focusing on model selection and uncertainty estimation in ecological applications, with an exttt{mgcv}-centric perspective.

Findings

Bayesian methods are integral to exttt{mgcv} and GAM theory.

Insights into model selection and uncertainty estimation in GAMs.

Practical guidance for ecologists using Bayesian and frequentist GAMs.

Abstract

Generalized additive models (GAMs) are a commonly used, flexible framework applied to many problems in statistical ecology. GAMs are often considered to be a purely frequentist framework (`generalized linear models with wiggly bits'), however links between frequentist and Bayesian approaches to these models were highlighted early on in the literature. Bayesian thinking underlies many parts of the implementation in the popular R package \texttt{mgcv} as well as in GAM theory more generally. This article aims to highlight useful links (and differences) between Bayesian and frequentist approaches to smoothing, and their practical applications in ecology (with an \texttt{mgcv}-centric viewpoint). Here I give some background for these results then move onto two important topics for quantitative ecologists: term/model selection and uncertainty estimation.

Figures

Equations

$$g({\mu }_i)={\mathbf{a}}_i^{⊺}\mathbf{\theta }+s_1(x_{1i})+s_2(x_{2i})+s_3(x_{3i},x_{4i})+\dots$$

$$s(x)=\sum _{k=1}^K{\beta }_k{b}_k(x),$$

$$\hat{\mathbf{\beta }}=\underset{\mathbf{\beta }}{\mathrm{argmax}}\left\{l(\mathbf{\beta })-\sum _{m=1}^M{\lambda }_m{\mathbf{\beta }}^{⊺}{\mathbf{S}}_m\mathbf{\beta }\right\},$$

$${\mathcal{L}}_p(\mathbf{\beta },\mathbf{\lambda })=\mathcal{L}(\mathbf{\beta })\exp \left(-{\mathbf{\beta }}^{⊺}{\mathbf{S}}_{\mathbf{\lambda }}\mathbf{\beta }\right).$$

$$\mathbb{E}\left({\text{CPUE}}_i\right)=\exp \left[{\beta }_0+s(x_i,y_i)+s({\text{Surface}}_i)+s({\text{Bottom}}_i)+s({\text{Depth}}_i)\right],$$

Full text

Centre for Research into Ecological & Environmental Modelling and School of Mathematics & Statistics, University of St Andrews, St Andrews, Fife, Scotland

1 Introduction

Flexible modelling of responses for a variety of distributions (binary, count, bounded, continuous) is an indispensable tool for quantitative ecologists. Generalized additive models [27, GAMs; e.g.,] provide one possible modelling framework in these situations. Common applications include species distribution modelling [8], abundance estimation [15], dose-response modelling [17], movement [1], ecosystem health [2] and more. In each case what is important is incorporating the structure of the data and/or data collection process into the model (be that the form of relationships, spatial correlation, blocking effects, etc). Informally this structure can be thought of as imposing some prior on how we would like the terms in the model to behave. In this article, I regurgitate some results from the statistical literature emphasizing this (Section 2) and then show how these tools can be used (or are already used) by those engaged in ecological modelling (Section 3).

Though the term “GAM” has significant baggage regarding the fitting method and type of terms, it really just describes the form of the linear predictor in the model (terms add together) and the response distribution. For example, a model may look like:

[TABLE]

where $\mu_{i}\equiv\mathbb{E}(Y_{i})$ and $Y_{i}\sim\text{EF}(\mu_{i},\phi)$ where $Y_{i}$ ($i=1,\ldots,n$) is the response and $\text{EF}(\mu_{i},\phi)$ indicates an exponential family distribution with mean $\mu_{i}$ and scale parameter $\phi$. $\bm{a}_{i}^{\intercal}$ is a vector of e.g., slopes and intercept terms, where $\bm{\theta}$ are their associated coefficients. The $s_{j}$ are “smooth” functions of one or more of the covariates ($x_{1i},x_{2i},x_{3i},x_{4i}$, etc). This definition can be adapted to Generalized Additive Mixed Models (GAMMs) and Generalized Linear Mixed Models (GLMMs), as we will see below.

The smooth terms are what makes GAMs interesting and useful, as compared to a generalized linear model. In a very general sense, they are constructed from sums of simple basis functions [6, e.g.,]: we make a complicated function by summing smaller, less complicated basis functions. In general for some smooth $s$ of covariate $x$ we have the following decomposition:

[TABLE]

where $b_{k}$ are fixed basis functions (with maximum complexity or basis dimension $K$) and $\beta_{k}$ are coefficients to be estimated. This basis function approach is extremely flexible, so to avoid overfitting we penalize the flexibility of each smooth term according to its wiggliness (so although $K$ may be set large, the fitted model has much smaller effective degrees of freedom, EDF). Generally such a penalty will be an integral of squared derivatives of $s$ (since derivatives measure the changes in $s$). The penalty can be written in the form $\sum_{m=1}^{M}\lambda_{m}\bm{\beta}^{\intercal}\bm{S}_{m}\bm{\beta}$ where $\bm{\beta}$ is a vector of coefficients, $\bm{S}_{m}$ is a matrix of the fixed parts of the penalty related to the $b_{k}$s and $\lambda_{m}$ are smoothing parameters to be estimated that control the influence of the penalty (where the $\bm{S}_{m}$s are padded with zeros so the sum forms a block matrix). Writing the penalty in this way means that we can compute the $\bm{S}_{m}$s once and during fitting the penalty is calculated by matrix multiplication only. Note that multiple penalty terms can correspond to a single smooth or multiple smooths may share a single smoothing parameter so $M$ is not necessarily the number of $s$ terms in the model.

We want to estimate model parameters that describe the data best, in the sense that we want to draw lines (or surfaces etc) that are close to the data but do not interpolate them. Formally we can set this up as a penalized log-likelihood [9]:

[TABLE]

where $l$ is the log-likelihood and there are $M$ smoothing parameters to estimate. Here were are trading-off between fit (high log-likelihood giving smooths close to the data) and penalty (large penalty indicating smooths are too wiggly).

Conditional on the $\lambda_{m}$s, estimation of $\hat{\bm{\beta}}$ in (3) is relatively straight-forward and the problem can be attacked with penalized iteratively re-weighted least squares (PIRLS) as for a GLM [27, section 6.1.1]. Estimating both $\hat{\bm{\beta}}$ and $\hat{\bm{\lambda}}$ is more complicated, as the smoothing parameters will constrain the values that $\hat{\bm{\beta}}$ can take. There have been various proposals for methods to fit such models in a frequentist framework and these fall into two categories: prediction error minimizing methods (like generalized cross-validation (GCV) or Akaike’s information criterion (AIC)) or likelihood based methods like (restricted maximum likelihood and marginal likelihood (REML/ML)) [26, see]. Prediction error minimizing methods have been shown to overfit (undersmooth) at finite sample sizes [20, for GCV;], especially when errors are correlated [10, for AIC;], so REML/ML have become the preferred methods. REML and ML cast the smooth functions as random effects [22] and smoothing parameters as variance parameters, so we can think of $\boldsymbol{\beta}$ as being a multivariate normal random effect with a variance proportional to $\boldsymbol{\lambda}$ (with structure imposed by the penalty matrix/matrices).

GAMs are often taught as an extension of the linear model: adding wiggles to make a (G)LM more flexible, often as a more principled step forward from adding polynomial terms. Smooth is often a synonym for spline, but there are many possible model terms that can be specified as basis functions subject to a sum of penalties: so-called “basis-penalty smoothers”. Taking covariates as-is (i.e., using the identity function) and having a diagonal penalty matrix implies a simple random effects model. Including more structure in the penalty, one can obtain more complex hierarchical random effects models. Having the penalty encode spatial information about a graphical structure gives (Gaussian) Markov random fields. We can also use tensor products of terms to construct multidimensional interaction-type effects, allowing for different units to be used for each covariate (i.e., anisotropy). Here I use the word ‘smooth’ to include all these possible flexible model terms.

When we talk about adding these smooth functions to our model, we tend to concentrate on equations like (1), looking at the mean effects of including smooths rather than thinking about the penalty. We usually view the penalty as a way of constraining our fit, stopping it from being too wiggly and ensuring that our model does not overfit. We can also think of the basis-penalty as the consequence of the problem definition, we have chosen them due to what we know about the dependencies and structures in the data (or data collection process). In a loose sense this is setting-up a prior on the kinds of functions we want to fit. The rest of this article investigates this idea further and explores some useful applications in ecology.

2 Bayesian interpretations

We can quickly get to a convenient Bayesian formulation by exponentiating the objective function in (3) [27, section 5.8], which gives us a penalized likelihood $\mathcal{L}_{p}$:

[TABLE]

We can then appeal to Bayes theorem and think of the exponential term as proportional to a multivariate normal distribution with mean zero and think of $\mathbf{S}_{\boldsymbol{\lambda}}$ (defined as $\mathbf{S}_{\lambda}=\sum_{m}\lambda_{m}\mathbf{S}_{m}$) as a prior precision matrix.

2.1 Specifying priors

By using a smooth term for a given covariate in our model we are specifying that observations which are close to each other in covariate space have similar values, the response varies smoothly (according to some measure of smoothness) and that the true function we seek to estimate is more likely to be smooth than wiggly. The Bayesian formulation allows us to be more explicit about these beliefs [27, section 5.8].

(4) says $\bm{\beta}\sim N(\bm{0},\bm{S}_{\lambda}^{-})$, where $\bm{S}_{\lambda}^{-}$ is the (pseudo)inverse of $\bm{S}_{\lambda}$. Large penalty entries in $\mathbf{S}$ correspond to wiggly basis functions (we want to penalize those more strongly) which, when inverted, give small variances (our prior is that basis function’s coefficient is close to zero): this makes explicit our belief that smoother functions are more likely than wiggly ones [30].

Since often some of the elements of $\boldsymbol{\beta}$ are not penalized (e.g., slope or intercept terms, which do not have derivatives), this leads to improper priors as there are no constraints on the value of the slope/intercept for those terms. In this case the pseudoinverse of $\bm{S}_{\lambda}$ is required. Some basis-penalty smoothers do lead to proper priors for all elements of $\boldsymbol{\beta}$ (e.g., the P-spline approach of [11]) and generally a proper prior can be found for any smooth by using the methods in Section 3.1. Identifiability constraints [27, section 5.4.1] that need to be imposed on the model (e.g., that there is only one intercept in the model) may also lead to proper priors [13].

Various different basis function-penalty combinations available in the literature express different priors on how we want our model terms to behave. For example, cyclic smoothers give us terms which ‘match’ up to a set number of derivatives at the start/end of the data and can be useful for temporal/seasonal effects. Many solutions have been proposed to the issue of spatial smoothers in areas with complex coastlines. [32] propose the soap film smoother, which estimates a boundary smooth which constrains the values inside the boundary. These kinds of effects are generated by transforms of the covariates (application of the basis functions) and the prior (penalty) on the corresponding coefficients.

REML and ML are referred to as empirical Bayes methods. When we take the random effects interpretation of $\boldsymbol{\beta}$, we can think of this as a prior and our fit criterion assesses the likelihood of the data given the implied prior on $\boldsymbol{\beta}$, as in (4). Here ‘empirical’ indicates that there is no prior for the smoothing parameters. Taking a fully Bayesian approach, it is common to put a vague gamma prior on each element of $\boldsymbol{\lambda}$ or uniform priors on their logarithm [31]. Specifying priors on smoothing/variance parameters can be tricky [24], this is especially the case for smoothing parameters as the true values of the smoothing parameter(s) could be infinite if e.g., the true smooth is in fact linear. It can also be hard to come-up with informed priors about smoothing parameters, as we often do not have a direct interpretation of their values.

As is often the case with richly-parameterized models, we can end up with many more coefficients than data. Though this is not necessarily a problem in a Bayesian setting, we may want to reduce the number of coefficients for computational/efficiency reasons. In a spline context, this tends to be done by limiting the number of basis functions ($K$ in (2)) to just those at knots: locations where the spline can change direction. The knot placement/number and basis complexity/dimension are usually linked such that there is one basis function per knot. Since $\bm{S}_{\lambda}^{-}$ involves basis functions , the number of basis functions or number of knots and their placement will affect the posterior. Effects of placement can be mitigated to some extent by over-specifying the number of knots/basis functions and allow wigglyness to be dictated by the smoothing parameter [19, 27, section 5.9]. Eigen-based approaches like thin-plate regression splines [29], make placement data-based in cases where regular grids are computationally taxing. Other related approaches include the use of triangulation-based techniques to optimize placement based on data locations [12].

2.2 Obtaining posteriors

For a fully Bayesian (FB) approach we formulate a likelihood and attach priors to the smoothing parameters $\boldsymbol{\lambda}$, as well as the model coefficients $\bm{\beta}$. We could then use MCMC to obtain a posterior. There are many software implementations which can achieve this, so here I only list R packages specifically tailored to GAMs: mgcv::jagam, which implements translation between mgcv and JAGS [31] or brms [4] which implements most mgcv models in Stan [5]. Dedicated software packages such as BayesX [3] can also be used. If one wishes to avoid MCMC, integrated nested Laplace approximations (INLA; e.g., via the R-INLA package) could be used [21, 28].

As discussed above, if we take an empirical Bayes (EB) view of the world and do not put priors on $\boldsymbol{\lambda}$, we can still obtain posteriors for $\bm{\beta}$, conditional on $\boldsymbol{\lambda}$. For computational efficiency the Laplace approximation is often used here; both R-INLA (used in ‘empirical Bayes mode’) and mgcv use this approach.

Using either EB or FB approaches we can get to the posterior marginal for $\bm{\beta}$: $\bm{\beta}|\bm{y},\bm{\lambda}\sim N(\hat{\bm{\beta}},\bm{V}_{\bm{\beta}})$ where for the Gaussian likelihood case $\bm{V}_{\bm{\beta}}=(\bm{X}^{\intercal}\bm{X}+\bm{S}_{\lambda})^{-1}\sigma^{2}$ and for the exponential family the expression is approximate and we have $\bm{V}_{\bm{\beta}}=(\bm{X}^{\intercal}\bm{W}\bm{X}+\bm{S}_{\lambda})^{-1}\phi$, where $\sigma^{2}$ is a variance parameter, $\phi$ is a scale parameter and $\mathbf{W}$ is a weight matrix [27, section 6.10]. For FB we can obtain a posterior for $\bm{\lambda}$ and an unconditional posterior for $\bm{\beta}$. [33] propose a correction to $\bm{V}_{\bm{\beta}}$ to account for the uncertainty in the smoothing parameter(s) using a Taylor expansion to approximate the extra variability when using EB.

In practice we can fit our models using EB methods (such as using REML/ML in mgcv) then sample from their posteriors. As noted, this is as straightforward as plugging the mean coefficient estimates and covariance matrix into a multivariate normal random number generator when using a Gaussian likelihood, though in the exponential family case one may have to use a Metropolis-Hasting sampler and propose from a t-distribution to get reasonable results (e.g., using mgcv::gam.mh).

3 Some examples

The Bayesian results above lead to some useful applications. Here I highlight a couple of the more commonly-used ones. To illustrate these techniques, data from the NOAA Alaska Fisheries Science Center’s groundfish assessment program (available at https://www.afsc.noaa.gov/RACE/groundfish/survey_data/default.htm) was used. The survey consists of summer bottom trawls at set of stations from 1982 through to 2018 and are shown in Figure Supplementary Figures. Response was catch per unit effort (CPUE; measured as individuals per hectare). Location (recorded as latitude/longitude but projected for analysis), date, surface temperature, bottom temperature (both recorded during the trawl, in degrees Celsius) and bathymetry (recorded in metres) were available as covariates. See [25] (and references therein) for further details of the survey. The examples below are not intended to be a serious analyses of the data. Data were downloaded from the NOAA AFSC website and processed for this analysis. The scripts used to process the data for this analysis are available at https://github.com/dill/RACE_bering_sea and code for each of the examples is available in the online Supplementary Material for this article and at https://github.com/dill/bayes-gam-explainer.

3.1 Term selection

We begin by fitting a model to the CPUE data for walleye pollock (Gadus chalcogrammus) in the Eastern Bering sea for 2010 only. The model includes a bivariate smooth of location and then univariate smooths of surface temperature, bottom temperature and bathymetry. Expected CPUE was modelled as

[TABLE]

where $i$ indexes the stations. CPUE was assumed to follow a Tweedie distribution with a log link [23, see e.g.,]. To model CPUE, we may not need all of the covariates: space ($x,y$), bottom depth (Depth), bottom temperature (Bottom) and surface temperature (Surface). Rather than using hypothesis testing for term selection, here I apply two different priors to remove terms during model fitting, following [13].

As described in section 2.1, the prior placed on $\bm{\beta}$ can be improper due to rank deficiency in $\mathbf{S}$. This means that there are e.g., linear or intercept terms that are not penalized. We refer to these terms as being in the nullspace of the penalty (the rest of the terms being referred to as the range space). We can make our priors proper by simply adding an extra penalty term to the model for the nullspace components of each term, equivalently we add an extra prior for these terms (the ‘double penalty’ approach of [13]). This is achieved by eigendecomposing the penalty matrix, $\mathbf{S}=\mathbf{U}\bm{\Lambda}\mathbf{U}^{\intercal}$. We can then form the additional penalty matrix $\mathbf{S}^{\ast}=\mathbf{U}^{\ast}\mathbf{U}^{\ast\intercal}$ where $\mathbf{U}^{*}$ is a matrix of eigenvectors corresponding to the zero entries on the diagonal of $\bm{\Lambda}$. Our original penalty $\mathbf{S}$ stays as-is, as the components in $\mathbf{S}^{*}$ do not have an effect (since their entries in $\boldsymbol{\Lambda}$ are (almost) zero). This approach is implemented as the select=TRUE option in mgcv::gam, and includes one additional smoothing parameter for each smooth term in the model, corresponding to each term’s nullspace. Alternatively, one can form a basis where the terms that lie in the nullspace have a shrinkage penalty applied to them by simply adding a small value to their corresponding diagonal entries of $\bm{\Lambda}$ so that the resulting penalty matrix is not rank-deficient (the ‘shrinkage’ approach of [13]; implemented as the cs and ts bases in mgcv). Equivalent to extending the prior on $\boldsymbol{\beta}$ to the terms in the nullspace.

These two approaches lead to rather different interpretations of how wigglyness should be penalized, or rather: the structure of the priors on the smooths. The shrinkage approach’s prior is that coefficients on wigglier terms should be reduced first, only once they have all been reduced to almost zero can the coefficients for the nullspace terms start to be reduced. So as the smoothing parameter increases (variance scaling decreases) the model goes from very wiggly, to just the terms in the nullspace (e.g., back to a linear model), to having no effect (estimated as zero). The double penalty approach prior makes no assumption about how much to smooth the nullspace components relative to the other parts of the smooth, the priors are separate (different smoothing parameters). This means that the nullspace components can be removed before the rest of the model, since there is a smoothing parameter for each part.

We can fit (5) in mgcv and see what the differences are between the results using these different prior specifications. Comparing the results from fitting all terms as thin plate regression splines (no selection), using the double penalty approach (with thin plate splines) and shrinkage (a modified thin plate spline) revealed that the two term selection methods completely removed the surface temperature term from the model. When no selection method was used, the surface temperature term remained as a linear term in the model (though it was not different from 0 according to an F-test). Figure 1 compares the resulting univariate smooth terms, though there are some minor differences the other smooths remain the same between the three models (though this is not guaranteed).

3.2 Uncertainty around smooth terms

From Section 2.2, we could use the posterior of $\bm{\beta}$ to generate possible parameters then use these to generate possible smooths. From these simulated smooths, we could then consider pointwise intervals over the range of the covariate to build percentile confidence (or credible) bands. Black lines in Figure 2 shows 1000 posterior samples of the smooth of depth for the shrinkage model in the previous section (black lines), and their 95% quantiles are the bounds of the blue band. We can take a shortcut and rather than simulating, we know that each smooth can be written as a linear combination ($s(x_{i})=\mathbf{X}_{i}\boldsymbol{\beta}$ for a model with a single smooth in it). We can then use a delta method-like approach to build point-wise credible intervals as $\hat{s}(x_{i})\pm z_{\alpha/2}\sqrt{v_{i}}$, where $\hat{s}$ is our estimated smooth, $v_{i}$ is the variance of the smooth at point $x_{i}$ and $z_{\alpha/2}$ is the usual appropriate value from a normal CDF (justification for these intervals was developed in [18] and expanded to the GAM case in [14]). These intervals have good frequentist across-the-function properties: that is, a 95% credible interval has close to 95% coverage, when averaged over the whole function. There may be over and under coverage at the peaks and troughs of the function as we know less about the exact turning points than we do about the function on the way to that turning point . The red band in Figure 2 shows these intervals, though is almost indistinguishable from the green band showing the quantile method. Since these intervals have good coverage and tell us about the whole function (by the across-the-function property), we can use them to test the hypothesis $H_{0}:s(x)=0\;\forall x$ —whether a term should be dropped from the model because it has no effect (the $p$-values presented in output of mgcv::summary).

We use the posterior samples in Figure 2 simply to calculate the blue band in the figure here but they can be useful beyond this. Simulating from the posterior of smooth terms (via simulation from the posterior of $\boldsymbol{\beta}$, conditional or not on $\boldsymbol{\lambda}$) can potentially reveal interesting properties of the fitted smooth which are not reflected in the plotted bands.

3.3 Posterior simulation/parametric bootstrap

Sometimes we want more than just uncertainty around individual terms in the model, we want to know about uncertainty either in the model’s predictions or summary statistics generated from predictions. Since we can simulate from the posterior of the model parameters, we can use those parameters to calculate functions of the simulated parameters. Calculating summary statistics on the results to obtain uncertainties about those quantities. This is particularly powerful as it allows us to calculate uncertainty about any function of the predictions (including transformations which are non-linear), avoiding potentially tricky derivations needed to obtain analytical expressions for the variance. A general algorithm [27, section 7.2.7] is as follows:

Let $B$ be the number of samples to generate. 2. 2.

Form $\mathbf{L}_{p}$, the matrix that maps the model covariates to the linear predictor (the prediction equivalent of the design matrix). 3. 3.

For $b$ in $1,\ldots,B$:

1. (a)

Simulate $\bm{\beta}_{b}$ from the posterior of $\hat{\bm{\beta}}$. 2. (b)

Calculate the linear predictor $\bm{\eta}_{b}=\mathbf{L}_{p}\bm{\beta}_{b}$. 3. (c)

Apply the inverse link function, $g$, so $\bm{y}_{b}=g^{-1}(\bm{\eta}_{b})$. 4. (d)

Calculate and store the required summary of $\bm{y}_{b}$. 4. 4.

Perform inference on the $B$ summaries (e.g., calculating empirical variance, percentile intervals, etc).

As an example of where we need to take summaries of non-linear functions of the linear predictor, we can fit a spatio-temporal model to all years (1982-2017) of walleye pollock CPUE data. Our model is then $\mathbb{E}\left(\text{CPUE}_{it}\right)=\exp\left[s(x_{it},y_{it},t)\right]$ (where $t$ indicates the year of the trawl at station $i$). Now the smooth $s()$ is constructed as a tensor product of a two dimensional thin plate spline smooth of $x$ and $y$, and a one dimensional thin plate spline smooth of $t$. The tensor product allows us to create multivariate interactions using simpler effects, in this case projecting our 2D spatial model through time to create a 3D spatio-temporal model [27, see e.g.,].

We want to obtain a time series of total predicted abundance at the stations per year. So we need to predict at each year and sum over the stations (i.e., space) at 3.(d) in the above algorithm. Figure 3 shows the predictions for the model, made by summing the predictions over space for each year, and their uncertainty. For simplicity, abundance is calculated by summing over the grid of all trawl station locations, it might be more appropriate to sum over a finer spatial grid and since surveys are in the summer only, abundance estimates were only made once per year (hence the piecewise linear nature of Figure 3). The green band in Figure 3 was generated by taking 10,000 samples with a Metropolis-Hastings sampler (as implemented in mgcv::gam.mh) with a burn-in of 5000 and thinning at every 10th sample. Increasing the spatial or temporal resolution involves modifying $\mathbf{L}_{p}$ and recalculating steps 3.(c), 3.(d) and 4. (simulation from the posterior does not have to be repeated).

4 Discussion

This article has highlighted the Bayesian interpretation of generalized additive models (specifically as implemented in mgcv), which are often thought of as a frequentist method. The article has emphasized that “GAM” only describes a model structure and that there are alternative ways to fit and interpret these models than are usually discussed. Work on Bayesian interpretations has shown that there are many ways which these links can exploited in practice for applied statistical work, they are not only of mathematical interest.

The Bayesian interpretations discussed here have been helpful to construct more reasonable estimates of uncertainty (including smoothing parameter uncertainty) and in order to understand how to construct confidence intervals that have good coverage properties. These links can surely be used further to develop other new methodology and enhance our understanding of the models that we fit. This approach has already been exploited to show that the stochastic partial differential equation approach proposed by [12] can be viewed as a basis-penalty smoother and implemented in mgcv [16]. Coming from the other direction, [7] expand on the idea of Bayesian regularisation and its interpretation, deriving corresponding priors for ridge regression, lasso, $L_{p}$ regularization, elastic net, etc. It is a shame that these conceptual links have not been better recognized and exploited further; even a very popular textbook [22] describes the mixed model representation of the GAM as a “convenient fiction”.

The jagam and brms packages allow ecologists to quickly build models using familiar syntax very similar to that for linear models, then transplant these into whatever fully Bayesian computation system they prefer. The models fitted in Section 3 could be fitted in e.g., JAGS or Nimble, using jagam to create necessary code. The main difference between those models and the ones presented here would be the priors on the smoothing parameters, which are not terribly interesting in these cases. Where these ideas really shine are in allowing smooths to be included as linear predictors for parameters in e.g., fully Bayesian occupancy or mark-recapture models. A general strategy that might be useful is using the GAM as a spatial distribution process for the study species, but building more complex observation processes (possibly from multiple data sources) in the native fully Bayesian package. In this way, the complex spatial structure is automatically generated and custom code is only required to stitch this part to the observation processes. Domain-specific knowledge can then be harnessed during prior specification, while allowing for extremely flexible model forms.

Moving beyond mere computational convenience and harnessing the broader Bayesian framework implicit in this modelling strategy can help increase understanding and synthesis, as well as providing further modelling extensions within a familiar framework.

Acknowledgements

This work was funded by OPNAV N45 and the SURTASS LFA Settlement Agreement, and being managed by the U.S. Navy’s Living Marine Resources program under Contract No. N39430-17-C-1982. Richard Glennie and Ian Durbach provided extremely useful feedback on an early draft.

Supplementary Figures

{sidewaysfigure}

Plot of the raw CPUE in space per year for walleye pollock. Empty grey circles indicate sampling effort but no catch, colours indicate catch per unit effort. Data from the NOAA Alaska Fisheries Science Center’s groundfish assessment program.