Understanding the stochastic partial differential equation approach to smoothing

01/21/2020 ∙ by David L Miller, et al. ∙ 0

Correlation and smoothness are terms used to describe a wide variety of random quantities. In time, space, and many other domains, they both imply the same idea: quantities that occur closer together are more similar than those further apart. Two popular statistical models that represent this idea are basis-penalty smoothers (Wood, 2017) and stochastic partial differential equations (SPDE) (Lindgren et al., 2011). In this paper, we discuss how the SPDE can be interpreted as a smoothing penalty and can be fitted using the R package mgcv, allowing practitioners with existing knowledge of smoothing penalties to better understand the implementation and theory behind the SPDE approach.



There are no comments yet.


page 12

page 19

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Data collected over space or time are often obtained with the desire to elicit an underlying pattern. The stochastic partial differential equation (SPDE) approach introduced by Lindgren et al. (2011) and implemented in the R-INLA software package (Rue et al., 2009) is a flexible, efficient method to analyse such data. Despite this, wider application is inhibited by two obstacles. First, the methods are presented using mathematical concepts and terms more usually found in applied mathematics and physics, making it difficult for practitioners in other fields to understand and adapt these methods to their own needs. Second, available software implementations are difficult to customise without high-level technical knowledge, limiting application to only those models available in the software or specially requested from software developers.

Here, we aim to mitigate these two issues: we describe () how the SPDE model can be interpreted as a basis-penalty smoother, a modelling framework more familiar to practitioners who use smoothing techniques (Wood, 2017), and () how software to fit these smoothers (e.g., mgcv), regularly extended and customised for application, can be used to fit SPDE models or, to go further, used to incorporate SPDE methods into larger models.

In this paper, we consider the following situation. Let

be a random variable observed at location

or time , depending on the domain. A statistical model for is constructed in three components or terms: . The first component, , is the fixed effect, often a linear combination of observed covariates with unknown parameters. The third component, , represents the measurement error or unstructured error, often for unknown parameter and every location . The second component is a stochastic process, representing the structured dependence among observations: observations made closer together in time or space are more likely to be similar than those further apart. A mathematically convenient and flexible model for this component is a Gaussian process (GP) with mean zero and covariance function . The covariance function quantifies how related two values of are at two locations. For fixed locations , the value of the GP at these locations, , are multivariate Gaussian with mean zero and covariance matrix with entry . We can extend this formulation to non-Gaussian responses by using a link function, , so the response is then modelled with a specified distribution and a mean .

An example of this kind of data might be a time series of counts. The left panel of Figure 1 show human cases of campylobacterosis (a common form of food poisoning, often originating in under-cooked poultry) in northern Québec, every 28 days from 1 January 1990 to 31 October 2000. We may expect a given time period’s count to be similar to its neighbours (e.g. due to seasonal variation), so our aim is to build a model that can capture this dependence. Using the above formulation, we model the counts as Poisson where represents time, is the number of cases at time , is a function of time representing the underlying process and is the appropriate inverse link function. Dependence structures become more complex when we move to a spatial domain. The right panel of Figure 1 shows remotely-sensed chlorophyll A levels in the Aral sea, derived from satellite data. In this case we expect that pixels close to each other have similar chlorophyll A levels. We now have represent a location in space, and the chlorophyll A level at that location, is modelled by , where now is a 2-dimensional stochastic process and . We revisit these examples in Section 4, below.

Figure 1: Examples of data with underlying dependence between observations. Left shows counts of campylobacterosis infections in northern Québec, summarized every 28 days from 1 January 1990 to 31 October 2000. Right shows the raw chlorophyll A in the Aral sea from the SeaWIFS satellite. In both cases we can build a model that takes into account the structure in the data.

Including in the model raises two issues: how to specify the covariance function and, once specified, how to fit the model. There are many possible solutions to these questions, including the SPDE and smoothing penalty approaches, and each uses different theoretical and numerical approximations; however, there is a common element: for observation locations , each method aims to define the covariance between these locations, by constructing an approximation to the precision matrix, defined as the inverse of the covariance matrix (). The precision matrix is a fundamental quantity required to fit these models (Simpson et al., 2012). The size of and makes the necessary computations expensive, in particular, if one of these matrices needs to be inverted so as to compute the other. It is in trying to avoid this computational burden that approximations are used.

The SPDE approach provides a method to approximate without the common substantial computational burden. The SPDE is an equation to be solved. Solutions to this equation are stochastic processes whose covariance structure is chosen to satisfy the relationship the SPDE specifies. The SPDE approach involves finding an SPDE whose solutions have the covariance structure, and implied precision matrix, that is desired for . Lindgren et al. (2011) show how to find an approximate solution to the SPDE by representing as a sum of basis functions multiplied by coefficients; this provides a computationally efficient way to compute : Lindgren et al. (2011) show that the coefficients of these basis functions form a Gaussian Markov random field, for which methods for fast computation of precision matrices already exist (Rue and Held, 2005). These computations make it possible to fit models quickly using integrated nested Laplace approximations (INLA; Rue et al., 2009).

Rue et al. (2017) is a comprehensive review of INLA that defines the class of latent Gaussian models with additive linear predictors which INLA is designed to fit. They also introduce Gaussian Markov random fields and their properties that lead to efficient computation, a key feature of the SPDE approach. Bakka et al. (2018) review the use of INLA for spatial modelling with a focus on the SPDE approach. They give an intuition for the method by introducing the notion of a “discretised” differential operator and describing the finite element methods that are used to solve the SPDE (Brenner and Scott, 2007; Bakka, 2018). See also Krainski et al. (2019) for a collection of worked examples of modelling with SPDEs using R-INLA, Wood (2019) for an approach to nested Laplace approximations without sparse Gaussian Markov random field structures, and Blangiardo and Cameletti (2015) for a comprehensive textbook on spatio-temporal modelling with R-INLA. We note that the R-INLA implementation of the SPDE approach has been applied in a wide variety of domains such as spatial epidemiology (Arab, 2015), species distribution mapping (de Rivera et al., 2018), spatial point processes (Simpson et al., 2016; Yuan et al., 2017; Soriano-Redondo et al., 2019) and environmental science (Huang et al., 2017), to name just a few examples. Our presentation here differs from the above resources in that we explicitly draw links with another well-known modelling framework.

The basis-penalty smoothing approach (Wood, 2017) is similar to the SPDE approach: the function is a sum of basis functions multiplied by coefficients. Rather than specify an SPDE and deduce a covariance structure between the coefficients, a smoothing penalty is used to induce correlation between the coefficients. This penalty measures how smooth is in its domain; intuitively, if changes more smoothly then values of at nearby locations are more correlated. Jointly optimising a measure of fit (sum of squares or log-likelihood) and smoothing penalty leads to an optimal curve, the smoothing spline (Wahba, 1990). This is a well-established approach with several excellent introductory resources (Hastie and Tibshirani, 1990; Ramsay and Silverman, 2005; Wood, 2017) and has been applied in many spatio-temporal modelling contexts (recent examples include Wood et al. (2017); Simpson (2018); Pedersen et al. (2019))

There is a direct correspondence between smoothing splines and stochastic processes (Kimeldorf and Wahba, 1970)

: the smoothing spline is a minimum variance unbiased linear estimator of the posterior mean of the stochastic process. For a stochastic process with a given covariance function, there is a corresponding SPDE and smoothing penalty such that one can estimate the posterior mean of

using the SPDE approach or the basis-penalty approach: both methods estimate the same quantity with the only differences being in numerical approximations and terminology. This means that the SPDE can be interpreted as a smoothing penalty and vice-versa.

This equivalence has been confirmed by Fahrmeir and Lang (2001), Lindgren and Rue (2008), and Yue et al. (2014) who show how basis-penalty smoothers in a Bayesian framework can be interpreted within the SPDE paradigm. Simpson et al. (2012) remark that the SPDE formulation is useful because it provides those with a background in physics or applied mathematics a way to understand and apply the model. In contrast, less emphasis has been placed on discussing this equivalence the other way around: SPDE methods can be formulated as basis-penalty smoothers. The SPDE formulation can seem opaque and fundamentally different for those unfamiliar with the mathematical concepts used. For this reason, showing the approach within the familiar smoothing framework demystifies the workings of the model and allows researchers in other fields to understand, adapt, and use the methods. We note that our approach is aligned with the general aim of emphasising links between Gaussian processes and the reproducing kernel Hilbert spaces theory that underpins the basis-penalty smoothing approach (Kanagawa et al., 2018), although here we take a more applied perspective.

Our aim in this paper is to show that the SPDE model as introduced by Lindgren et al. (2011) (usually fitted using R-INLA) can be described as a basis-penalty smoother and fitted using mgcv. To do this, we first describe the SPDE method for those unfamiliar with the mathematical concepts used, highlighting the key steps in the method. Afterward, we show the equivalences and differences between the SPDE method and the analogous basis-penalty smoother.

2 The SPDE approach

2.1 What is an SPDE?

A stochastic partial differential equation involves stochastic processes and differential operators. Examples of differential operators () are the first derivative, the second derivative, the gradient operator in two-dimensions or the Laplacian in two dimensions. Combinations of these are also differential operators, e.g., such that for a function . Here, we consider only linear differential operators, that is, is a linear combination of derivatives of , of different orders. Differential operators of stochastic processes can be treated similarly to those applied to ordinary functions, there is one key difference that we will highlight below. Overall, an SPDE states that the differential of a function

is equal to some known stochastic process, most commonly the white noise process,

. The white noise process is completely uncorrelated and is a Normal random variable with mean zero and finite variance for every .

In general, the SPDE states that for some differential operator . A stochastic process, , is called a solution to the SPDE if it satisfies this equation. Consider an example, let be the first derivative of the function. The SPDE therefore states that the first derivative of has mean zero and finite variance at every point; furthermore, it states that the value of the derivative at points and are uncorrelated for all . Approximately, this means that for a small and point , where is a Normal random variable. Consider if the SPDE has a parameter such that such that controls the variance in the white noise process. This means that changes in are more variable when is reduced and less variable for higher . In other words, the parameters of the SPDE control the smoothness of . It is important to note that here the term “smoothness” is not used in a mathematical way, meaning differentiability, nor in a strictly statistical way, referring to correlation range, but in a qualitative way—when we speak of differentiability or correlation we shall use those terms explicitly.

For a given , the mathematical form of the solution to the SPDE is known: where is a function you can derive given you know . The function is called Green’s function; in the appendix (Proposition ) we show how this function is derived from . Intuitively, acts as a weighting function such that the value of the stochastic process at is a weighted sum over the white noise process; this is called a convolution. Suppose were set to give infinite weight to distance , , and zero elsewhere, for , then : is just the white noise process, completely uncorrelated. Alternatively, if gave equal weight to all distances, e.g., for all , then would be constant, perfectly correlated. Between these two extremes are weighting functions that reproduce correlations over different ranges. It can be shown that the covariance function is given by , see appendix (Proposition ) for the derivation.

In summary, the solutions to the SPDE have a covariance structure that is induced by the choice of . This means that one could describe a system using an SPDE and then deduce the associated covariance function from it. The power of the SPDE approach is realised by doing the opposite: find a that induces the covariance function that you want. The power of finding the SPDE corresponding to a desired covariance function is that the precision matrix can be efficiently computed using the SPDE.

2.2 Solving the SPDE

The SPDE involves applying a differential operator to a stochastic process, , but this cannot be done in the same way as when you apply to a known function. This is because is random and, in many cases, realisations of will not be suitably differentiable. For example, the Brownian motion stochastic process has a derivative equal to the white noise process, but it is also known that simulated trajectories of Brownian motion are nowhere differentiable. is a convenient shorthand way to think about the SPDE, but technically, the SPDE only has meaning when stated in an integral form. That is, means that we require for every function with compact support. The function is often called the test function. For brevity, let and so the integral form is . The notation is called the inner product of , it has many nice mathematical properties, including being linear, that is for functions and constants .

In the integral form, the equation makes sense because any stochastic process can be integrated, but not every one can be differentiated. By requiring the equation to hold for every , we require the left-hand stochastic process and the right-hand process to have the same integral, no matter how we average over space. For example, if the stochastic processes were one-dimensional, we could split the real line into intervals and select a function to be one on this interval and zero outside. Since the integral equation must hold for all such functions, we therefore require to have the same average value as on each and every interval.

Given an SPDE, Lindgren et al. (2011) show how to derive an approximate solution using the finite element method. The domain (e.g., time or space) is split into “elements”, e.g., a grid or a triangulation, often called a mesh. To each point on this mesh, a basis function is associated. The solution to the SPDE is then a weighted sum of the basis functions and random variables : .

The integral form of the SPDE then implies that for any function , . We cannot, however, check this equation for infinitely many test functions , so instead we restrict to only testing with the functions that can be written in our chosen basis. As is a linear operator, this is equivalent to solving the system of equations for every . This system can be written as a matrix equation: where has entry and has entry .

To summarise, the SPDE is written in an integral form, sometimes using inner products, since stochastic processes are well defined when integrated but not when differentiated. Given this, the solution is represented in a chosen basis. The integral form is then solved by considering only test functions within that basis. This leads to a matrix equation involving the coefficients , the matrix

, and the random vector

. The random vector has known distribution, because it depends only on the basis functions and the white noise process:

has a multivariate Gaussian distribution with mean zero and a precision matrix

where has entry . It follows from that where . The SPDE is therefore a way to specify a prior for .

This provides an approximate solution to the SPDE. For example, given an SPDE, one can use the finite element method to compute and therefore simulate from a multivariate Gaussian distribution with precision . The function would then be a realisation from a stochastic process which is a solution to the SPDE, a stochastic process with the covariance structure implied by .

2.3 Matérn SPDE

The focus of Lindgren et al. (2011) and the covariance function most commonly used in the R-INLA software is the Matérn covariance function. The Matérn covariance function is considered a flexible model for the dependencies found in real world observations: it has the form where are parameters, is the modified Bessel function of the second kind, and is the dimension of the domain. Figure 2 shows realisations from two stochastic processes with Matérn covariance functions in one-dimension, one with a longer correlation range than the other.

It is difficult to fit models with this covariance structure due to the computational issues mentioned above. Lindgren et al. (2011) apply the finite element method to approximate stochastic processes with Matérn covariance (a comparison of the notation used in Section 2.2 and that used in Lindgren et al. (2011) is given in the appendix, Section 5). To do this, they present the differential operator that corresponds to this covariance function: where .

When , this is called a fractional differential operator; for this paper, we consider only the case when and so is again a linear differential operator. In practice, is poorly identified and difficult to estimate from data, so its value is often assumed to be fixed (Zhang, 2004).

Lindgren et al. (2011) solve the SPDE using the finite element method. By deriving the weighting function and computing the covariance from this SPDE, Whittle (1954) shows that the solutions have Matérn covariance, as desired. In other words, the precision matrix computed from the finite element method is an approximation to the precision matrix one would obtain if you computed the variance-covariance matrix with the Matérn covariance function and then, at great computational cost, inverted this matrix. Figure 2 shows a subview of the approximate precision matrix: the matrix is mostly filled with zeroes, with non-zero values occurring on three bands down the diagonal. This is an example of a sparse matrix, computations with these matrices are efficient because many of the computations ordinarily required can be omitted as it is known the matrix is mostly zeroes.

To use the finite element method one must chose a mesh, a grid or triangulation, over the domain and a basis to define on this mesh. The default choice in R-INLA is to use a regular grid in 1D (or a constrained Delaunay triangulation in 2D) to produce a mesh and then define piecewise linear basis functions (specifically, linear B-spline basis functions) on this mesh.

Figure 2: Two functions, one smooth (long-range correlation, dashed line, open circle data) and one rough (short-range correlation, solid line, filled circle data) (top plot), their Matérn correlation functions (bottom left plot, same line types) and the first 11 rows and columns of an example approximate Gaussian Markov Random field precision matrix (bottom right plot, darker shade indicates higher absolute value, each row and column corresponds to a data point location).

3 The basis-penalty approach

3.1 What is a basis-penalty smoother?

The basis-penalty approach refers to models where is assumed to have the form for basis functions and parameters . A model is then assumed for the observations given this form for to provide a measure of fit, the log-likelihood . Alternatively, the sum-of-squares can be used as a measure of fit. Optimising to obtain leads to a function that, given

is large enough, will interpolate the observed data: capturing the noise in the observations as well as the underlying signal (such overfitting stops us from making inference on the signal). A smoothing penalty,

, is subtracted from the log-likelihood to penalize functions that are too wiggly. The smoothing parameter, , controls the extent of the penalization (a larger value of leads to a smoother ).

The estimates for are defined to be those that optimise the joint measure of fit and smoothness, the penalised likelihood: . This involves estimating both the optimal smoothing parameter and coefficients . In practice, REstricted Maximum Likelihood (REML; Wood, 2011) is used to to do this.

There are several choices for the smoothing penalty. Most are defined using a differential operator . For example, in one dimension, the smoothing penalty , i.e., where is the second derivative, is often used. For this penalty, functions with rapidly changing gradients are penalised while functions with constant gradient, straight lines, have no penalty. In higher dimensions, the thin-plate spline (Wood, 2003) is often used with penalty: for two-dimensions. This penalty takes the total variation in the gradient of including the interaction between the coordinates. The penalty for smoothing splines takes the form for some chosen differential operator (see Yue and Speckman (2010) and Yue et al. (2014) who show this for the thin plate spline penalty). This can also be written as an inner product .

When , the penalty based on the differential operator can be written in matrix form: where is a matrix with entry .

In summary, a basis-penalty smoother is specified by selecting a basis, e.g., a B-spline basis of specified order, and a smoothing penalty. The parameters are then estimated by optimising the penalised likelihood: .

3.2 Connection between SPDE and penalty

Rewriting the penalized log-likelihood as a likelihood we obtain . A Bayesian interpretation of the penalised likelihood as proportional to a posterior implies that is an improper prior for (Silverman, 1985; Wood, 2017). Since

is proportional to a multivariate normal distribution with mean zero and precision matrix

, the penalized likelihood is equivalent to assigning the prior .

The connection between the SPDE approach and the basis-penalty approach can now be made clear. It can be shown that for a given differential operator , the approximate precision matrix for the SPDE is the same as the precision matrix computed using the smoothing penalty (appendix, Proposition 3).

This connection has two implications. First, it means that the differences between the basis-penalty approach and the SPDE finite element approximation, when using the same basis and differential operator, are differences in implementation only, as both should lead to the same approximate precision matrix. Second, the connection means that any SPDE of the form can be understood and interpreted as a smoothing penalty of the form , and vice-versa.

3.3 Matérn penalty

The SPDE specified in Lindgren et al. (2011) has the differential operator . Given the connection described above, this can be interpreted as a smoothing penalty: . This penalty is different from those considered above because it contains two smoothing parameters: and . This offers it more flexibility. The penalty can still, however, be interpreted as such: it is a trade-off between the value of the function and the second derivative in each direction. As is increased, the value of increases, meaning that can be higher, the function be less smooth, while keeping the penalty the same. Alternatively, can be described as the inverse correlation range: higher values of lead to less smooth functions meaning values of the function become less correlated. The smoothing parameter controls the overall smoothness of .

The Matérn penalty can be written in matrix form as above, but for computational convenience, it is first broken into three parts: . Notice that it appears that the Laplacian has been replaced with the gradient operator : this relationship holds here using Green’s first identity and the Neumann boundary condition, see Bakka et al. (2018) for more detail. This leads to the smoothing matrix where are all matrices with entries , , and , respectively. All of these matrices are sparse and so computation of the smoothing penalty, , is computationally efficient. The matrix is equal to the matrix computed using the finite element method (Appendix, Proposition 3).

3.4 Fitting the Matérn SPDE in mgcv

The mgcv R package allows the specification of new basis-penalty smoothers by writing new “smooth.construct” functions which build an appropriate design matrix (containing evaluations of the basis functions), penalty matrices and other optional components. Within this framework we can fit the SPDE model in mgcv providing a smooth.construct.spde.smooth.spec constructor. R-INLA provides helper functions to construct the required design and penalty matrices. Here we sketch an algorithm for setting-up SPDE models in mgcv.

Given we have a response and covariates in an matrix , we construct our model as follows.

  1. Create a mesh using INLA::inla.mesh.1d or INLA::inla.mesh.2d.

  2. Calculate , and using INLA::inla.mesh.fem (c1, g1 and g2, respectively).

  3. We need to connect the basis representation of to the observation locations. At present contains the value of at each mesh point, not at each observation location. A matrix multiplication is used to project the values at all mesh points to the observations locations, it is called the projection matrix (found using INLA::inla.spde.make.A). The full design matrix is then given by combining the fixed effects design matrix and the contribution for , .

  4. Having constructed the design matrix and penalty matrices, use REML to find optimal , and subject to the penalty matrix . (Parametrisation for this model in mgcv is given in Supplementary Material section 4.)

As REML is an empirical Bayes procedure, we expect point estimates for to coincide for the procedure outlined above and R-INLA. A uniform prior is implied for the smoothing parameters ( or ); R-INLA

allows for similar estimation by just using the modes of the hyperparameters

and (the int.strategy="eb" option). Proper priors could be used if step (4), above, was replaced by an MCMC scheme.

4 Examples

We now compare the SPDE and basis-penalty models applied to three example datasets. We fitted the SPDE Matérn model in both R-INLA and mgcv. Code for these examples is available as supplementary material.

4.1 Campylobacterosis cases in Québec

Ferland et al. (2006) analyse a time series of (human) cases of campylobacterosis in northern Québec, with observations every 28 days from 1 January 1990 to 31 October 2000 (140 observations). We modelled the number of infections as a function of time, using a Poisson response and a link function. The model is fitted using three approaches () a Matérn basis-penalty smoother with 50 degree 2 B-splines, fitted with mgcv; () a Matérn SPDE for with a finite element basis of 50 degree 2 B-splines and penalized complexity priors (Simpson et al., 2017) on smoothing parameters, fitted with R-INLA; () a basis-penalty smoother with penalty equal to the integral of the squared second derivative, using 50 degree 2 B-splines fitted using mgcv.

Figure 3: Campylobacterosis cases modelled using: a Matérn basis-penalty smoother fitted with mgcv (top), a Matérn SPDE fitted with R-INLA (middle), a B-spline basis-penalty smoother fitted using mgcv (bottom).

Fitted models are shown in Figure 3. Results from the R-INLA and mgcv SPDE implementations are very similar. This is supported by the similarity in the estimated hyperparameters ( and for R-INLA, and and for mgcv). The squared second derivative penalty B-spline fit from mgcv is smoother than those from the SPDE-based methods.

4.2 Aral sea chlorophyll

Moving to a 2-dimensional smoothing problem, we consider remotely sensed () chlorophyll from the Aral sea collected by the NASA SeaWifs satellite over a series of 8 day observation periods. The 496 observations used here are averages (from 1998 to 2002) of the 38th observation period. Data were taken from the gamair package (dataset aral) and consist of spatial coordinates and logarithm of chlorophyll concentration.

We built a mesh using fmesher::meshbuilder (Supplementary Figure 2) and generated two-dimensional degree 1 B-splines. We consider the model for location with no fixed effects. We fitted the Matérn model using the SPDE and penalty approaches in R-INLA and mgcv. For the R-INLA model, penalized complexity priors were used.

There was little visual difference with good agreement in the predictions (Supplementary Figure 3 shows the posterior mode and percentile credible surfaces for these models and Supplementary Figure 4 gives a graphical comparison). Hyperparameter estimates were similar: and for R-INLA and and for mgcv

. To investigate differences between the two models we took (1000) samples from the posterior of each model and looked at the differences between pairs of realisations on a per-cell basis. Plots of the mean of these differences and their standard deviations are shown in Figure

4. The mean plot shows structure to the differences in the models, though differences are relatively small (range of chlorophyll A values in original data: 1.905–19.275). This is to be expected if the models produce similar predicted values to each other, which are consistent through each realisation.

Figure 4: Chlorophyll in the Aral sea example. Left shows mean difference in predictions and right shows standard deviation of the difference in predictions between SPDE models fitted using mgcv and R-INLA.

4.3 MODIS land surface temperatures

To compare the two approaches on a large data set, we now turn to land surface temperature data collected by the Terra instrument on the MODIS satellite. The data consist of a grid of measurements in the latitude range 34.29519–37.06811 and longitude range -95.91153—91.28381 on August 4, 2016. The training data (105,569 observations) as defined in Heaton et al. (2018) were used to fit the model, but a significantly simpler mesh was used (see Supplementary Figure 5). Following from Heaton et al. (2018) we assumed a Gaussian response for temperature and fitted a 2-dimensional SPDE model on latitude and longitude.

As the dataset was quite large, we used the bam (“big additive model”) function in mgcv to fit the SPDE model (additionally discretizing covariate values for efficient storage and computation; Wood et al., 2017). The SPDE fitted by R-INLA used the empirical Bayes integration strategy. We timed the fitting for both approaches (ignoring mesh setup, which was shared across methods), taking only the time for inla and bam to run. The mgcv model was slightly faster (4.71 minutes versus 5.23 in R-INLA running on a Windows server using 1 core of a Xeon Gold 6152 at 2.1GHz with 512Gb RAM). Supplementary Figure 6 compares the predictions from the two methods, the largest absolute difference between predictions is 4.761, which is small compared to the range of the data.

5 Discussion

We have drawn links between two approaches to fitting the same model: the stochastic partial differential equation method as implemented in R-INLA and the basis-penalty smoothing approach as fitted in mgcv. This paper aims to make accessible what is equivalent between the approaches, what is a matter of choice, and what is fundamentally different. Yue et al. (2014) show how splines can be specified using the SPDE approach, benefitting those familiar with SPDEs. Here, we do the opposite for the benefit of those familiar with the (penalized likelihood/empirical Bayes) GAM framework. Supplementary Figure 1 is a flow diagram showing the parallels between the smoothness and correlation approaches we have discussed.

Similarities between many smoothing techniques can be drawn. Smoothing splines, kriging, Gaussian Markov random fields, and SPDEs approximate similar models, but their explanations make it difficult for practitioners to appreciate their commonalities and determine precisely what is a necessary and what is a coincidental association. Taking the precision matrix as the common currency between these methods, a modelling framework emerges:

  1. Choose a covariance model: explicitly, as in kriging, through the smoothness penalty as with smoothing splines, or with an SPDE;

  2. Approximate the precision matrix : reduce dimension (fixed rank kriging, thin plate splines) or induce sparsity in (B-splines, SPDE);

  3. Draw approximate inference using a software implementation: e.g., with mgcv, MCMC (e.g., Stan (Carpenter et al., 2017); JAGS, (Plummer, 2017)), R-INLA (Rue et al., 2009), lme4 (Bates et al., 2015) or TMB (Kristensen et al., 2016).

This paper is an example of comparing two methods according to this framework. Doing so for other smoothing methods will allow alternative modelling approaches to be compared on the grounds of genuine differences: in the covariance function, in the approximation for , in the estimation procedure, or, simply, in the software implementation.


The authors wish to thank the editor and two anonymous reviewers for their constructive comments on the first submission of the manuscript. They also thank Finn Lindgren for extremely helpful input and Simon Wood for the suggestion of the smoothing parameter parameterisation. The authors also thank Steve Buckland, David Borchers and Joe Watson for their comments on the manuscript.


  • Arab (2015) Arab, A. (2015). Spatial and spatio-temporal models for modeling epidemiological data with excess zeros. International Journal of Environmental Research and Public Health 12(9), 10536–10548.
  • Bakka (2018) Bakka, H. (2018). How to solve the stochastic partial differential equation that gives a Matérn random field using the finite element method. arXiv preprint arXiv:1803.03765.
  • Bakka et al. (2018) Bakka, H., H. Rue, G.-A. Fuglstad, A. Riebler, D. Bolin, J. Illian, E. Krainski, D. Simpson, and F. Lindgren (2018). Spatial modeling with r-inla: A review. Wiley Interdisciplinary Reviews: Computational Statistics 10(6), e1443.
  • Bates et al. (2015) Bates, D., M. Mächler, B. Bolker, and S. Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software 67(1).
  • Blangiardo and Cameletti (2015) Blangiardo, M. and M. Cameletti (2015). Spatial and Spatio-temporal Bayesian Models with R-INLA. Wiley, Chichester, UK.
  • Brenner and Scott (2007) Brenner, S. and R. Scott (2007). The mathematical theory of finite element methods, Volume 15. Springer, New York, NY, USA.
  • Carpenter et al. (2017) Carpenter, B., A. Gelman, M. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. Brubaker, J. Guo, P. Li, and A. Riddell (2017). Stan: A Probabilistic Programming Language. Journal of Statistical Software 76(1), 1–32.
  • de Rivera et al. (2018) de Rivera, O. R., M. Blangiardo, A. López-Quílez, and I. Martín-Sanz (2018). Species distribution modelling through Bayesian hierarchical approach. Theoretical Ecology, 1–11.
  • Fahrmeir and Lang (2001) Fahrmeir, L. and S. Lang (2001). Bayesian inference for generalized additive mixed models based on Markov random field priors. Journal of the Royal Statistical Society: Series C (Applied Statistics) 50(2), 201–220.
  • Ferland et al. (2006) Ferland, R., A. Latour, and D. Oraichi (2006, November). Integer-Valued GARCH Process. Journal of Time Series Analysis 27(6), 923–942.
  • Hastie and Tibshirani (1990) Hastie, T. and R. Tibshirani (1990). Generalized Additive Models.

    Number 43 in Monographs on Statistics and Applied Probability. Chapman and Hall, London, UK.

  • Heaton et al. (2018) Heaton, M. J., A. Datta, A. O. Finley, R. Furrer, J. Guinness, R. Guhaniyogi, F. Gerber, R. B. Gramacy, D. Hammerling, M. Katzfuss, F. Lindgren, D. W. Nychka, F. Sun, and A. Zammit-Mangion (2018, December). A case study competition among methods for analyzing large spatial data. Journal of Agricultural, Biological and Environmental Statistics.
  • Huang et al. (2017) Huang, J., B. P. Malone, B. Minasny, A. B. McBratney, and J. Triantafilis (2017). Evaluating a Bayesian modelling approach (INLA-SPDE) for environmental mapping. Science of the Total Environment 609, 621–632.
  • Kanagawa et al. (2018) Kanagawa, M., P. Hennig, D. Sejdinovic, and B. K. Sriperumbudur (2018, July). Gaussian Processes and Kernel Methods: A Review on Connections and Equivalences. arXiv preprint arXiv:1807.02582. arXiv: 1807.02582.
  • Kimeldorf and Wahba (1970) Kimeldorf, G. S. and G. Wahba (1970). A correspondence between Bayesian estimation on stochastic processes and smoothing by splines. The Annals of Mathematical Statistics 41(2), 495–502.
  • Krainski et al. (2019) Krainski, E., V. Gómez-Rubio, H. Bakka, A. Lenzi, D. Castro-Camilo, D. Simpson, F. Lindgren, and H. Rue (2019). Advanced Spatial Modeling with Stochastic Partial Differential Equations Using R and INLA. CRC Press, Boca Raton, FL, USA.
  • Kristensen et al. (2016) Kristensen, K., A. Nielsen, C. W. Berg, H. Skaug, and B. M. Bell (2016). TMB: Automatic differentiation and Laplace approximation. Journal of Statistical Software 70(5), 1–21.
  • Lindgren and Rue (2008) Lindgren, F. and H. Rue (2008). On the second-order random walk model for irregular locations. Scandinavian journal of statistics 35(4), 691–700.
  • Lindgren et al. (2011) Lindgren, F., H. Rue, and J. Lindström (2011). An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73(4), 423–498.
  • Pedersen et al. (2019) Pedersen, E. J., D. L. Miller, G. L. Simpson, and N. Ross (2019). Hierarchical generalized additive models in ecology: an introduction with mgcv. PeerJ 7, e6876.
  • Plummer (2017) Plummer, M. (2017). JAGS Version 4.3.0 user manual.
  • Ramsay and Silverman (2005) Ramsay, J. and B. W. Silverman (2005). Functional Data Analysis (2 ed.). Springer Series in Statistics. New York: Springer-Verlag.
  • Rue and Held (2005) Rue, H. and L. Held (2005). Gaussian Markov Random Fields: Theory and Applications, Volume 104 of Monographs on Statistics and Applied Probability. Chapman & Hall, London, UK.
  • Rue et al. (2009) Rue, H., S. Martino, and N. Chopin (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 71(2), 319–392.
  • Rue et al. (2017) Rue, H., A. Riebler, S. H. Sørbye, J. B. Illian, D. P. Simpson, and F. K. Lindgren (2017). Bayesian computing with inla: A review. Annual Review of Statistics and Its Application 4(1), 395–421.
  • Silverman (1985) Silverman, B. W. (1985). Some aspects of the spline smoothing approach to non-parametric regression curve fitting. Journal of the Royal Statistical Society. Series B (Methodological) 47(1), 1–52.
  • Simpson et al. (2016) Simpson, D., J. B. Illian, F. Lindgren, S. H. Sørbye, and H. Rue (2016). Going off grid: computationally efficient inference for log-Gaussian Cox processes. Biometrika 103(1), 49–70.
  • Simpson et al. (2012) Simpson, D., F. Lindgren, and H. Rue (2012). In order to make spatial statistics computationally feasible, we need to forget about the covariance function. Environmetrics 23(1), 65–74.
  • Simpson et al. (2017) Simpson, D., H. Rue, A. Riebler, T. G. Martins, and S. H. Sørbye (2017). Penalising model component complexity: a principled, practical approach to constructing priors. Statistical Science 32(1), 1–28.
  • Simpson (2018) Simpson, G. L. (2018). Modelling Palaeoecological Time Series Using Generalised Additive Models. Frontiers in Ecology and Evolution 6.
  • Soriano-Redondo et al. (2019) Soriano-Redondo, A., C. M. Jones-Todd, S. Bearhop, G. M. Hilton, L. Lock, A. Stanbury, S. C. Votier, and J. B. Illian (2019). Understanding species distribution in dynamic populations: a new approach using spatio-temporal point process models. Ecography 42(6), 1092–1102.
  • Wahba (1990) Wahba, G. (1990). Spline models for observational data. SIAM, Philadephia, PA, USA.
  • Whittle (1954) Whittle, P. (1954). On stationary processes in the plane. Biometrika 41(3-4), 434–449.
  • Wood (2003) Wood, S. N. (2003). Thin plate regression splines. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 65(1), 95–114.
  • Wood (2011) Wood, S. N. (2011). Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73(1), 3–36.
  • Wood (2017) Wood, S. N. (2017). Generalized Additive Models. An Introduction with R (2nd ed.). Texts in Statistical Science. CRC Press, Boca Raton, FL, USA.
  • Wood (2019) Wood, S. N. (2019). Simplified integrated nested Laplace approximation. Biometrika.
  • Wood et al. (2017) Wood, S. N., Z. Li, G. Shaddick, and N. H. Augustin (2017, July). Generalized Additive Models for Gigadata: Modeling the U.K. Black Smoke Network Daily Data. Journal of the American Statistical Association 112(519), 1199–1210.
  • Yuan et al. (2017) Yuan, Y., F. E. Bachl, F. Lindgren, D. L. Borchers, J. B. Illian, S. T. Buckland, H. Rue, T. Gerrodette, et al. (2017). Point process models for spatio-temporal distance sampling data from a large-scale survey of blue whales. The Annals of Applied Statistics 11(4), 2270–2297.
  • Yue and Speckman (2010) Yue, Y. and P. L. Speckman (2010). Nonstationary spatial gaussian markov random fields. Journal of Computational and Graphical Statistics 19(1), 96–116.
  • Yue et al. (2014) Yue, Y. R., D. Simpson, F. Lindgren, and H. Rue (2014). Bayesian adaptive smoothing splines using stochastic differential equations. Bayesian Analysis 9(2), 397–424.
  • Zhang (2004) Zhang, H. (2004). Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics. Journal of the American Statistical Association 99(465), 250–261.