In this vignette we will present the `R-INLA`

implementation of
the rational SPDE approach. For theoretical details we refer the reader
to the Rational approximation with the
`rSPDE`

package vignette and to Xiong,
Simas, and Bolin (2022).

We begin by providing a step-by-step illustration on how to use our implementation. To this end we will consider a real world data set that consists of precipitation measurements from the Paraná region in Brazil.

After the initial model fitting, we will show how to change some parameters of the model. In the end, we will also provide an example in which we have replicates.

It is important to mention that one can improve the performance by
using the PARDISO solver. Please, go to https://www.pardiso-project.org/r-inla/#license to apply
for a license. Also, use `inla.pardiso()`

for instructions on
how to enable the PARDISO sparse library.

To illustrate our implementation of `rSPDE`

in `R-INLA`

we will consider a
dataset available in `R-INLA`

. This data has
also been used to illustrate the SPDE approach, see for instance the
book Advanced
Spatial Modeling with Stochastic Partial Differential Equations Using R
and INLA and also the vignette Spatial
Statistics using R-INLA and Gaussian Markov random fields. See also
Lindgren, Rue, and Lindström (2011) for
theoretical details on the standard SPDE approach.

The data consist of precipitation measurements from the Paraná region in Brazil and were provided by the Brazilian National Water Agency. The data were collected at 616 gauge stations in Paraná state, south of Brazil, for each day in 2011.

We will follow the vignette Spatial
Statistics using R-INLA and Gaussian Markov random fields. As
precipitation data are always positive, we will assume it is Gamma
distributed. `R-INLA`

uses the following parameterization of the Gamma distribution, \[\Gamma(\mu, \phi): \pi (y) =
\frac{1}{\Gamma(\phi)} \left(\frac{\phi}{\mu}\right)^{\phi} y^{\phi - 1}
\exp\left(-\frac{\phi y}{\mu}\right) .\] In this
parameterization, the distribution has expected value \(E(x) = \mu\) and variance \(V(x) = \mu^2/(\phi)\), where\(1/\phi\) is a dispersion parameter.

In this example \(\mu\) will be modeled using a stochastic model that includes both covariates and spatial structure, resulting in the latent Gaussian model for the precipitation measurements \[\begin{align} y_i\mid \mu(s_i), \theta &\sim \Gamma(\mu(s_i),c\phi)\\ \log (\mu(s)) &= \eta(s) = \sum_k f_k(c_k(s))+u(s)\\ \theta &\sim \pi(\theta) \end{align},\]

where \(y_i\) denotes the
measurement taken at location \(s_i\),
\(c_k(s)\) are covariates, \(u(s)\) is a mean-zero Gaussian Matérn
field, and \(\theta\) is a vector
containing all parameters of the model, including smoothness of the
field. That is, by using the `rSPDE`

model we will also be
able to estimate the smoothness of the latent field.

We will be using `R-INLA`

. To install `R-INLA`

go to R-INLA Project.

We begin by loading some libraries we need to get the data and build the plots.

```
library(ggplot2)
library(INLA)
library(splancs)
library(viridis)
```

Let us load the data and the border of the region

```
data(PRprec)
data(PRborder)
```

The data frame contains daily measurements at 616 stations for the year 2011, as well as coordinates and altitude information for the measurement stations. We will not analyze the full spatio-temporal data set, but instead look at the total precipitation in January, which we calculate as

`<- rowMeans(PRprec[, 3 + 1:31]) Y `

In the next snippet of code, we extract the coordinates and altitudes and remove the locations with missing values.

```
<- !is.na(Y)
ind <- Y[ind]
Y <- as.matrix(PRprec[ind, 1:2])
coords <- PRprec$Altitude[ind] alt
```

Let us build plot the precipitation observations using
`ggplot`

:

```
ggplot() +
geom_point(aes(
x = coords[, 1], y = coords[, 2],
colour = Y
size = 2, alpha = 1) +
), geom_path(aes(x = PRborder[, 1], y = PRborder[, 2])) +
geom_path(aes(x = PRborder[1034:1078, 1], y = PRborder[
1034:1078,
2
colour = "red") +
]), scale_color_viridis()
```

The red line in the figure shows the coast line, and we expect the distance to the coast to be a good covariate for precipitation. This covariate is not available, so let us calculate it for each observation location:

```
<- apply(spDists(coords, PRborder[1034:1078, ],
seaDist longlat = TRUE
1, min) ),
```

Now, let us plot the precipitation as a function of the possible covariates:

```
par(mfrow = c(2, 2))
plot(coords[, 1], Y, cex = 0.5, xlab = "Longitude")
plot(coords[, 2], Y, cex = 0.5, xlab = "Latitude")
plot(seaDist, Y, cex = 0.5, xlab = "Distance to sea")
plot(alt, Y, cex = 0.5, xlab = "Altitude")
```

`par(mfrow = c(1, 1))`

To use the `R-INLA`

implementation of the `rSPDE`

model we need to load the
functions:

`library(rSPDE)`

The `rSPDE`

-`INLA`

implementation is very
reminiscent of `R-INLA`

,
so its usage should be straightforward for `R-INLA`

users. For
instance, to create a `rSPDE`

model, one would use
`rspde.matern()`

in place of
`inla.spde2.matern()`

. To create an index, one should use
`rspde.make.index()`

in place of
`inla.spde.make.index()`

. To create the `A`

matrix, one should use `rspde.make.A()`

in place of
`inla.spde.make.A()`

, and so on.

The main differences when comparing the arguments between the
`rSPDE`

-`INLA`

implementation and the standard
SPDE implementation in `R-INLA`

, are the
`nu`

and `rspde.order`

arguments, which are
present in `rSPDE`

-`INLA`

implementation. We will
see below how use these arguments.

We can use `R-INLA`

for creating the mesh. Let us create a mesh which is based on a
non-convex hull to avoid adding many small triangles outside the domain
of interest:

```
<- inla.nonconvex.hull(coords, -0.03, -0.05, resolution = c(100, 100))
prdomain <- inla.mesh.2d(boundary = prdomain, max.edge = c(0.45, 1), cutoff = 0.2)
prmesh plot(prmesh, asp = 1, main = "")
lines(PRborder, col = 3)
points(coords[, 1], coords[, 2], pch = 19, cex = 0.5, col = "red")
```

We now create the \(A\) matrix, that
connects the mesh to the observation locations and then create the
`rSPDE`

model.

For this task, as we mentioned earlier, we need to use an
`rSPDE`

specific function, whose name is very reminiscent to
`R-INLA`

’s standard SPDE
approach, namely `rspde.make.A()`

(in place of `R-INLA`

’s
`inla.spde.make.A()`

). The reason for the need of this
specific function is that the size of the \(A\) matrix depends on the order of the
rational approximation. The details can be found in the introduction of
the Rational approximation with the
`rSPDE`

package vignette.

The default order is 2 for our covariance-based rational
approximation. As mentioned in the introduction of the Rational approximation with the `rSPDE`

package vignette, an approximation of order 2 in the
covariance-based rational approximation has approximately the same
computational cost as the operator-based rational approximation of order
1.

Recall that the latent process \(u\)
is a solution of \[(\kappa^2
I-\Delta)^{\alpha/2}(\tau u) = \mathcal{W},\] where \(\alpha = \nu + d/2\). We want to estimate
all three parameters \(\tau,\kappa\)
and \(\nu\), which is the default
option of

the `rSPDE`

-`INLA`

implementation. However, there
is also an option to fix the smoothness parameter \(\nu\) to some predefined value and only
estimate \(\tau\) and \(\kappa\). This will be discussed later.

In this first example we will assume we want a rational approximation
of order 2. To this end we can use the `rspde.make.A()`

function. Since we will assume order 2 and that we want to estimate
smoothness, which are the default options in this function, the required
parameters are simply the mesh and the locations:

`<- rspde.make.A(mesh = prmesh, loc = coords) Abar `

To set up an `rSPDE`

model, all we need is the mesh. By
default it will assume that we want to estimate the smoothness parameter
\(\nu\) and to do a covariance-based
rational approximation of order 2.

Later in this vignette we will also see other options for setting up
`rSPDE`

models such as keeping the smoothness parameter fixed
and/or increasing the order of the covariance-based rational
approximation.

Therefore, to set up a model all we have to do is use the
`rspde.matern()`

function:

`<- rspde.matern(mesh = prmesh) rspde_model `

Note that this function is very reminiscent of `R-INLA`

’s
`inla.spde2.matern()`

function. This is a pattern we have
tried to keep consistent in the package: All the `rSPDE`

versions of some `R-INLA`

function will
either replace `inla`

or `inla.spde`

or
`inla.spde2`

by `rspde`

.

`inla.stack`

Since the covariates are already evaluated at the observation
locations, we only want to apply the \(A\) matrix to the spatial effect and not
the fixed effects. We can use the `inla.stack()`

function.

The difference, however, is that we need to use the function
`rspde.make.index()`

(in place of the standard
`inla.spde.make.index()`

) to create the index.

If one is using the default options, that is, to estimate the
smoothness parameter \(\nu\) and to do
a rational approximation of order 2, the usage of
`rspde.make.index()`

is identical to the usage of
`inla.spde.make.index()`

:

`<- rspde.make.index(name = "field", mesh = prmesh) mesh.index `

We can then create the stack in a standard manner:

```
<- inla.stack(
stk.dat data = list(y = Y), A = list(Abar, 1), tag = "est",
effects = list(
c(
mesh.index,list(Intercept = 1)
),list(
seaDist = inla.group(seaDist)
)
) )
```

Here the observation matrix \(A\) is applied to the spatial effect and the intercept while an identity observation matrix, denoted by \(1\), is applied to the covariates. This means the covariates are unaffected by the observation matrix.

The observation matrices in \(A=list(Abar,1)\) are used to link the
corresponding elements in the effects-list to the observations. Thus in
our model the latent spatial field `mesh.index`

and the
intercept are linked to the log-expectation of the observations,
i.e. \(\eta(s)\), through the \(A\)-matrix. The covariates, on the other
hand, are linked directly to \(\eta(s)\). The `stk.dat`

object
defined above implies the following principal linkage between model
components and observations \[\eta(s) \sim A
x(s) + A \text{ Intercept} + \text{seaDist}.\] \(\eta(s)\) will then be used in the
observation-likelihood, \[y_i\mid
\eta(s_i),\theta \sim \Gamma(\exp(\eta (s_i)), c\phi).\]

We will build a model using the distance to the sea \(x_i\) as a covariate through an improper
CAR(1) model with \(\beta_{ij}=1(i\sim
j)\), which `R-INLA`

calls a random
walk of order 1.

```
<- y ~ -1 + Intercept + f(seaDist, model = "rw1") +
f.s f(field, model = rspde_model)
```

Here `-1`

is added to remove R’s implicit intercept, which
is replaced by the explicit `+Intercept`

from when we created
the stack.

To fit the model we proceed as in the standard SPDE approach and we
simply call `inla()`

.

```
<- inla(f.s,
rspde_fit family = "Gamma", data = inla.stack.data(stk.dat),
verbose = FALSE,
control.inla = list(int.strategy = "eb"),
control.predictor = list(A = inla.stack.A(stk.dat), compute = TRUE)
)
```

We can look at some summaries of the posterior distributions for the parameters, for example the fixed effects (i.e. the intercept) and the hyper-parameters (i.e. dispersion in the gamma likelihood, the precision of the RW1, and the parameters of the spatial field):

`summary(rspde_fit)`

```
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 3.73, Running = 11.7, Post = 0.0677, Total = 15.5
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 0.647 0.018 0.612 0.647 0.683 0.647 0
##
## Random effects:
## Name Model
## seaDist RW1 model
## field CGeneric
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision parameter for the Gamma observations 13.199 0.901 11.484
## Precision for seaDist 10587.150 8024.021 2455.436
## Theta1 for field -1.468 0.160 -1.788
## Theta2 for field 0.011 0.343 -0.638
## Theta3 for field -0.454 0.919 -2.202
## 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 13.180 15.031 13.162
## Precision for seaDist 8394.282 31882.128 5491.261
## Theta1 for field -1.467 -1.157 -1.460
## Theta2 for field 0.001 0.714 -0.038
## Theta3 for field -0.478 1.423 -0.571
##
## Marginal log-Likelihood: -1263.53
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
```

Let \(\theta_1 = \textrm{Theta1}\), \(\theta_2=\textrm{Theta2}\) and \(\theta_3=\textrm{Theta3}\). In terms of the SPDE \[(\kappa^2 I - \Delta)^{\alpha/2}(\tau u) = \mathcal{W},\] where \(\alpha = \nu + d/2\), we have that \[\tau = \exp(\theta_1),\quad \kappa = \exp(\theta_2), \] and by default \[\nu = 4\Big(\frac{\exp(\theta_3)}{1+\exp(\theta_3)}\Big).\] The number 4 comes from the upper bound for \(\nu\), which will be discussed later in this vignette.

In general, we have \[\nu = \nu_{UB}\Big(\frac{\exp(\theta_3)}{1+\exp(\theta_3)}\Big),\] where \(\nu_{UB}\) is the value of the upper bound for the smoothness parameter \(\nu\).

Another choice for prior for \(\nu\) is a truncated lognormal distribution and will also be discussed later in this vignette.

`rSPDE`

-`INLA`

resultsWe can obtain outputs with respect to parameters in the original
scale by using the function `rspde.result()`

:

```
<- rspde.result(rspde_fit, "field", rspde_model)
result_fit summary(result_fit)
```

```
## mean sd 0.025quant 0.5quant 0.975quant mode
## std.dev 0.233254 0.037119 0.167832 0.230773 0.313302 0.226317
## range 1.072870 0.385643 0.531836 0.999776 2.027850 0.873352
## nu 1.210560 0.571970 0.303711 1.145120 2.406620 0.888165
```

To create plots of the posterior marginal densities, we can use the
`gg_df()`

function, which creates
`ggplot2`

-friendly data frames. The following figure shows
the posterior marginal densities of the three parameters using the
`gg_df()`

function.

```
<- gg_df(result_fit)
posterior_df_fit
ggplot(posterior_df_fit) + geom_line(aes(x = x, y = y)) +
facet_wrap(~parameter, scales = "free") + labs(y = "Density")
```

This function is reminiscent to the `inla.spde.result()`

function with the main difference that it has the `summary()`

and `plot()`

methods implemented.

Let us now obtain predictions (i.e. do kriging) of the expected precipitation on a dense grid in the region.

We begin by creating the grid in which we want to do the predictions.
To this end, we can use the `rspde.mesh.projector()`

function. This function has the same arguments as the function
`inla.mesh.projector()`

, with the only difference being that
the `rSPDE`

version also has an argument `nu`

and
an argument `rspde.order`

. Thus, we proceed in the same
fashion as we would in `R-INLA`

’s standard SPDE
implementation:

```
<- c(150, 100)
nxy <- rspde.mesh.projector(prmesh,
projgrid xlim = range(PRborder[, 1]),
ylim = range(PRborder[, 2]), dims = nxy
)
```

This lattice contains 150 × 100 locations. One can easily change the
resolution of the kriging prediction by changing `nxy`

. Let
us find the cells that are outside the region of interest so that we do
not plot the estimates there.

`<- inout(projgrid$lattice$loc, cbind(PRborder[, 1], PRborder[, 2])) xy.in `

Let us plot the locations that we will do prediction:

```
<- projgrid$lattice$loc[xy.in, ]
coord.prd plot(coord.prd, type = "p", cex = 0.1)
lines(PRborder)
points(coords[, 1], coords[, 2], pch = 19, cex = 0.5, col = "red")
```

Now, there are a few ways we could calculate the kriging prediction. The simplest way is to evaluate the mean of all individual random effects in the linear predictor and then to calculate the exponential of their sum (since \(\mu(s)=\exp(\eta(s))\) ). A more accurate way is to calculate the prediction jointly with the estimation, which unfortunately is quite computationally expensive if we do prediction on a fine grid. However, in this illustration, we proceed with this option to show how one can do it.

To this end, first, link the prediction coordinates to the mesh nodes through an \(A\) matrix

`<- projgrid$proj$A[xy.in, ] A.prd `

Since we are using distance to the sea as a covariate, we also have to calculate this covariate for the prediction locations.

```
<- apply(spDists(coord.prd,
seaDist.prd 1034:1078, ],
PRborder[longlat = TRUE
1, min) ),
```

We now make a stack for the prediction locations. We have no data at
the prediction locations, so we set `y= NA`

. We then join
this stack with the estimation stack.

```
<- list(
ef.prd c(mesh.index, list(Intercept = 1)),
list(
long = inla.group(coord.prd[
,1
lat = inla.group(coord.prd[, 2]),
]), seaDist = inla.group(seaDist.prd)
)
)<- inla.stack(
stk.prd data = list(y = NA),
A = list(A.prd, 1), tag = "prd",
effects = ef.prd
)<- inla.stack(stk.dat, stk.prd) stk.all
```

Doing the joint estimation takes a while, and we therefore turn off
the computation of certain things that we are not interested in, such as
the marginals for the random effect. We will also use a simplified
integration strategy (actually only using the posterior mode of the
hyper-parameters) through the command
`control.inla = list(int.strategy = "eb")`

, i.e. empirical
Bayes.

```
<- inla(f.s,
rspde_fitprd family = "Gamma",
data = inla.stack.data(stk.all),
control.predictor = list(
A = inla.stack.A(stk.all),
compute = TRUE, link = 1
),control.compute = list(
return.marginals = FALSE,
return.marginals.predictor = FALSE
),control.inla = list(int.strategy = "eb")
)
```

We then extract the indices to the prediction nodes and then extract the mean and the standard deviation of the response:

```
<- inla.stack.index(stk.all, "prd")$data
id.prd <- rspde_fitprd$summary.fitted.values$mean[id.prd]
m.prd <- rspde_fitprd$summary.fitted.values$sd[id.prd] sd.prd
```

Finally, we plot the results:

```
# Plot the predictions
<- data.frame(x1 = coord.prd[,1],
pred_df x2 = coord.prd[,2],
mean = m.prd,
sd = sd.prd)
ggplot(pred_df, aes(x = x1, y = x2, fill = mean)) +
geom_raster() +
scale_fill_viridis()
```

Then, the std. deviations:

```
ggplot(pred_df, aes(x = x1, y = x2, fill = sd)) +
geom_raster() + scale_fill_viridis()
```

For this example we will simulate a data with replicates. We will use
the same example considered in the Rational
approximation with the `rSPDE`

package vignette (the only
difference is the way the data is organized). We also refer the reader
to this vignette for a description of the function
`matern.operators()`

, along with its methods (for instance,
the `simulate()`

method).

Let us consider a simple Gaussian linear model with 30 independent replicates of a latent spatial field \(x(\mathbf{s})\), observed at the same \(m\) locations, \(\{\mathbf{s}_1 , \ldots , \mathbf{s}_m \}\), for each replicate. For each \(i = 1,\ldots,m,\) we have

\[\begin{align} y_i &= x_1(\mathbf{s}_i)+\varepsilon_i,\\ \vdots &= \vdots\\ y_{i+29m} &= x_{30}(\mathbf{s}_i) + \varepsilon_{i+29m}, \end{align}\]

where \(\varepsilon_1,\ldots,\varepsilon_{30m}\) are iid normally distributed with mean 0 and standard deviation 0.1.

We use the basis function representation of \(x(\cdot)\) to define the \(A\) matrix linking the point locations to
the mesh. We also need to account for the fact that we have 30
replicates at the same locations. To this end, the \(A\) matrix we need can be generated by
`inla.spde.make.A()`

function. The reason being that we are
sampling \(x(\cdot)\) directly and not
the latent vector described in the introduction of the Rational approximation with the `rSPDE`

package vignette.

We begin by creating the mesh:

```
<- 200
m <- matrix(runif(m * 2), m, 2)
loc_2d_mesh <- inla.mesh.2d(
mesh_2d loc = loc_2d_mesh,
cutoff = 0.05,
offset = c(0.1, 0.4),
max.edge = c(0.05, 0.5)
)plot(mesh_2d, main = "")
points(loc_2d_mesh[, 1], loc_2d_mesh[, 2])
```

We then compute the \(A\) matrix, which is needed for simulation, and connects the observation locations to the mesh:

```
<- 30
n.rep <- inla.spde.make.A(
A mesh = mesh_2d,
loc = loc_2d_mesh,
index = rep(1:m, times = n.rep),
repl = rep(1:n.rep, each = m)
)
```

Notice that for the simulated data, we should use the \(A\) matrix from
`inla.spde.make.A()`

function.

We will now simulate a latent process with standard deviation \(\sigma=1\) and range \(0.1\). We will use \(\nu=0.5\) so that the model has an
exponential covariance function. To this end we create a model object
with the `matern.operators()`

function:

```
<- 0.5
nu <- 1
sigma <- 0.1
range <- sqrt(8 * nu) / range
kappa <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) * (4 * pi) * gamma(nu + 1)))
tau <- 2
d <- matern.operators(
operator_information mesh = mesh_2d,
nu = nu,
kappa = kappa,
sigma = sigma,
m = 2
)
```

More details on this function can be found at the Rational approximation with the rSPDE package vignette.

To simulate the latent process all we need to do is to use the
`simulate()`

method on the `operator_information`

object. We then obtain the simulated data \(y\) by connecting with the \(A\) matrix and adding the gaussian
noise.

```
set.seed(1)
<- simulate(operator_information, nsim = n.rep)
u <- as.vector(A %*% as.vector(u)) +
y rnorm(m * n.rep) * 0.1
```

The first replicate of the simulated random field as well as the observation locations are shown in the following figure.

```
<- inla.mesh.projector(mesh_2d, dims = c(100, 100))
proj
<- data.frame(x = proj$lattice$loc[,1],
df_field y = proj$lattice$loc[,2],
field = as.vector(inla.mesh.project(proj,
field = as.vector(u[, 1]))),
type = "field")
<- data.frame(x = loc_2d_mesh[, 1],
df_loc y = loc_2d_mesh[, 2],
field = y[1:m],
type = "locations")
<- rbind(df_field, df_loc)
df_plot
ggplot(df_plot) + aes(x = x, y = y, fill = field) +
facet_wrap(~type) + xlim(0,1) + ylim(0,1) +
geom_raster(data = df_field) +
geom_point(data = df_loc, aes(colour = field),
show.legend = FALSE) +
scale_fill_viridis() + scale_colour_viridis()
```

`## Warning: Removed 7599 rows containing missing values (`geom_raster()`).`

Let us then use the rational SPDE approach to fit the data.

We begin by creating the \(A\)
matrix and index with replicates, and the `inla.stack`

object. It is important to notice that since we have replicates we
should provide the `index`

and `repl`

arguments
for `rspde.make.A()`

function, and also the argument
`n.repl`

in `rspde.make.index()`

function. They
behave identically as in their `R-INLA`

’s counterparts,
namely, `inla.spde.make.A()`

and
`inla.make.index()`

.

```
<- rspde.make.A(
Abar.rep mesh = mesh_2d, loc = loc_2d_mesh, index = rep(1:m, times = n.rep),
repl = rep(1:n.rep, each = m)
)<- rspde.make.index(
mesh.index.rep name = "field", mesh = mesh_2d,
n.repl = n.rep
)
<- inla.stack(
st.dat.rep data = list(y = y),
A = Abar.rep,
effects = mesh.index.rep
)
```

We now create the model object.

`<- rspde.matern(mesh = mesh_2d, parameterization = "spde") rspde_model.rep `

Finally, we create the formula and fit. It is extremely important not
to forget the `replicate`

argument when building the formula
as `inla()`

function will not produce warning and might fit
some meaningless model.

```
<-
f.rep ~ -1 + f(field,
y model = rspde_model.rep,
replicate = field.repl
)<-
rspde_fit.rep inla(f.rep,
data = inla.stack.data(st.dat.rep),
family = "gaussian",
control.predictor =
list(A = inla.stack.A(st.dat.rep))
)
```

We can get the summary:

`summary(rspde_fit.rep)`

```
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 3.28, Running = 116, Post = 1.37, Total = 120
## Random effects:
## Name Model
## field CGeneric
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 91.48 3.263 85.76 91.25
## Theta1 for field -3.01 0.074 -3.19 -3.01
## Theta2 for field 3.10 0.015 3.07 3.10
## Theta3 for field -1.28 0.036 -1.37 -1.28
## 0.975quant mode
## Precision for the Gaussian observations 98.61 90.39
## Theta1 for field -2.83 -3.01
## Theta2 for field 3.14 3.10
## Theta3 for field -1.19 -1.28
##
## Marginal log-Likelihood: -4400.88
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
```

and the summary in the user’s scale:

```
<- rspde.result(rspde_fit.rep, "field", rspde_model.rep)
result_fit_rep summary(result_fit_rep)
```

```
## mean sd 0.025quant 0.5quant 0.975quant mode
## tau 0.0493901 0.00365111 0.0425381 0.0492587 0.056900 0.0490386
## kappa 22.2631000 0.32873300 21.6283000 22.2594000 22.919100 22.2519000
## nu 0.6512080 0.01848920 0.6156890 0.6509310 0.688295 0.6503330
```

```
<- data.frame(
result_df parameter = c("tau", "kappa", "nu"),
true = c(tau, kappa, nu),
mean = c(
$summary.tau$mean,
result_fit_rep$summary.kappa$mean,
result_fit_rep$summary.nu$mean
result_fit_rep
),mode = c(
$summary.tau$mode,
result_fit_rep$summary.kappa$mode,
result_fit_rep$summary.nu$mode
result_fit_rep
)
)print(result_df)
```

```
## parameter true mean mode
## 1 tau 0.08920621 0.04939006 0.04903857
## 2 kappa 20.00000000 22.26311581 22.25194969
## 3 nu 0.50000000 0.65120777 0.65033268
```

It is also possible to consider models in which \(\sigma\) (std. deviation) and \(\rho\) (range parameter) are non-stationary. One can also use the parameterization in terms of the SPDE parameters \(\kappa\) and \(\tau\).

An example of such a model is given in the vignette inlabru implementation of the rational SPDE approach.

`rSPDE`

-`INLA`

implementationWe will now discuss some of the arguments that were introduced in our
`R-INLA`

implementation
of the rational approximation that are not present in `R-INLA`

’s standard SPDE
implementation.

In each case we will provide an illustrative example.

When we fit a `rspde.matern()`

model we need to provide an
upper bound for the smoothness parameter \(\nu\). The reason for that is that the
sparsity of the precision matrix should be kept fixed during `R-INLA`

’s estimation and
the higher the value of \(\nu\) the
denser the precision matrix gets.

This means that the higher the value of \(\nu\), the higher the computational cost to fit the model. Therefore, ideally, want to choose an upper bound for \(\nu\) as small as possible.

To change the value of the upper bound for the smoothness parameter,
we must change the argument `nu.upper.bound`

. The default
value for `nu.upper.bound`

is 4. Other common choices for
`nu.upper.bound`

are 2 and 1.

It is clear from the discussion above that the smaller the value of
`nu.upper.bound`

the faster the estimation procedure will
be.

However, if we choose a value of `nu.upper.bound`

which is
too low, the “correct” value of \(\nu\)
might not belong to the interval \((0,\nu_{UB})\), where \(\nu_{UB}\) is the value of
`nu.upper.bound`

. Hence, one might be forced to increase
`nu.upper.bound`

and estimate again, which, obviously will
increase the computational cost as we will need to do more than one
estimation.

Let us illustrate by considering the same model we considered above
for the precipitation in Paraná region in Brazil and consider
`nu.upper.bound`

equal to 2, which is generally a good choice
for `nu.upper.bound`

.

We simply use the function `rspde.matern()`

with the
argument `nu.upper.bound`

set to 2:

`<- rspde.matern(mesh = prmesh, nu.upper.bound = 2) rspde_model_2 `

Since we are considering the default `rspde.order`

, the
\(A\) matrix and the mesh index objects
are the same as the previous ones. Let us then update the formula and
fit the model:

```
.2 <- y ~ -1 + Intercept + f(seaDist, model = "rw1") +
f.sf(field, model = rspde_model_2)
<- inla(f.s.2,
rspde_fit_2 family = "Gamma", data = inla.stack.data(stk.dat),
verbose = FALSE,
control.inla = list(int.strategy = "eb"),
control.predictor = list(A = inla.stack.A(stk.dat), compute = TRUE)
)
```

Let us see the summary of the fit:

`summary(rspde_fit_2)`

```
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 3.2, Running = 4.37, Post = 0.056, Total = 7.62
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 0.648 0.02 0.608 0.648 0.688 0.648 0
##
## Random effects:
## Name Model
## seaDist RW1 model
## field CGeneric
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision parameter for the Gamma observations 13.190 0.899 11.496
## Precision for seaDist 10874.412 8322.223 3003.789
## Theta1 for field -1.364 0.153 -1.643
## Theta2 for field 0.209 0.423 -0.554
## Theta3 for field -0.939 0.721 -2.478
## 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 13.165 15.036 13.122
## Precision for seaDist 8533.088 33204.409 5680.041
## Theta1 for field -1.371 -1.040 -1.402
## Theta2 for field 0.186 1.113 0.086
## Theta3 for field -0.908 0.382 -0.752
##
## Marginal log-Likelihood: -1262.51
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
```

Let us compare with the cost from the previous fit, with the default
value of `nu.upper.bound`

of 4:

```
# nu.upper.bound = 4
$cpu.used rspde_fit
```

```
## Pre Running Post Total
## 3.72995400 11.66820002 0.06769991 15.46585393
```

```
# nu.upper.bound = 2
$cpu.used rspde_fit_2
```

```
## Pre Running Post Total
## 3.19765210 4.37015700 0.05600905 7.62381816
```

We can see that the fit for `nu.upper.bound`

equal to 2
was considerably faster.

Finally, let us get the result results for the field and see the estimate of \(\nu\):

```
<- rspde.result(rspde_fit_2, "field", rspde_model_2)
result_fit_2 summary(result_fit_2)
```

```
## mean sd 0.025quant 0.5quant 0.975quant mode
## std.dev 0.258640 0.0403282 0.194084 0.253303 0.351538 0.241050
## range 1.350060 0.6313410 0.580263 1.196080 2.997810 0.953979
## nu 0.602964 0.2702980 0.159606 0.581452 1.175140 0.511953
```

To change the order of the rational approximation all we have to do
is set the argument `rspde.order`

to the desired value. The
current available possibilities are `1,2,3`

,…, up to
`8`

.

The higher the order of the rational approximation, the more accurate the results will be, however, the higher the computational cost will be.

The default `rspde.order`

of 2 is generally a good choice
and reasonably accurate. See the vignette Rational approximation with the rSPDE package
for further details on the order of the rational approximation and some
comparison with the Matérn covariance.

Let us fit the above model with the covariance-based rational
approximation of order `3`

. Since we are changing the order
of the rational approximation, that is, we are changing the
`rspde.order`

argument, we need to recompute the \(A\) matrix and the mesh index. Therefore,
we proceed as follows:

- We build a new model:

```
<- rspde.matern(mesh = prmesh,
rspde_model_order_3 rspde.order = 3,
nu.upper.bound = 2
)
```

- We create a new \(A\) matrix:

`<- rspde.make.A(mesh = prmesh, loc = coords, rspde.order = 3) Abar_3 `

- We create a new index:

```
.3 <- rspde.make.index(
mesh.indexname = "field", mesh = prmesh,
rspde.order = 3
)
```

Now the remaining steps are the same as before:

```
.3 <- inla.stack(
stk.datdata = list(y = Y), A = list(Abar_3, 1), tag = "est",
effects = list(
c(
.3,
mesh.indexlist(Intercept = 1)
),list(
long = inla.group(coords[, 1]),
lat = inla.group(coords[, 2]),
seaDist = inla.group(seaDist)
)
)
)
.3 <- y ~ -1 + Intercept + f(seaDist, model = "rw1") +
f.sf(field, model = rspde_model_order_3)
<- inla(f.s.3,
rspde_fit_order_3 family = "Gamma", data = inla.stack.data(stk.dat.3),
verbose = FALSE,
control.inla = list(int.strategy = "eb"),
control.predictor = list(A = inla.stack.A(stk.dat.3), compute = TRUE)
)
```

Let us see the summary:

`summary(rspde_fit_order_3)`

```
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 3.17, Running = 12.6, Post = 0.0637, Total = 15.8
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 0.486 0.014 0.459 0.486 0.513 0.486 0
##
## Random effects:
## Name Model
## seaDist RW1 model
## field CGeneric
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision parameter for the Gamma observations 13.213 0.900 11.505
## Precision for seaDist 11300.152 9917.588 2538.668
## Theta1 for field -1.435 0.162 -1.751
## Theta2 for field 0.080 0.385 -0.628
## Theta3 for field -0.308 0.951 -2.248
## 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 13.192 15.047 13.167
## Precision for seaDist 8410.016 37588.024 5232.242
## Theta1 for field -1.437 -1.111 -1.443
## Theta2 for field 0.062 0.885 -0.008
## Theta3 for field -0.283 1.496 -0.184
##
## Marginal log-Likelihood: -1263.38
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
```

We can see in the above summary that the computational cost
significantly increased. Let us compare the cost of having
`rspde.order=3`

and `nu.upper.bound=2`

with the
cost of having `rspde.order=2`

and
`nu.upper.bound=4`

:

```
# nu.upper.bound = 4
$cpu.used rspde_fit
```

```
## Pre Running Post Total
## 3.72995400 11.66820002 0.06769991 15.46585393
```

```
# nu.upper.bound = 2
$cpu.used rspde_fit_order_3
```

```
## Pre Running Post Total
## 3.17286706 12.60617900 0.06371093 15.84275699
```

One can check the order of the rational approximation by using the
`rational.order()`

function. It also allows another way to
change the order of the rational order, by using the corresponding
`rational.order<-()`

function.

The `rational.order()`

and
`rational.order<-()`

functions can be applied to the
`inla.rspde`

object, to the `A`

matrix and to the
`index`

objects.

Here to check the models:

`rational.order(rspde_model)`

`## [1] 2`

`rational.order(rspde_model_order_3)`

`## [1] 3`

Here to check the `A`

matrices:

`rational.order(Abar)`

`## [1] 2`

`rational.order(Abar_3)`

`## [1] 3`

Here to check the indexes:

`rational.order(mesh.index)`

`## [1] 2`

`rational.order(mesh.index.3)`

`## [1] 3`

Let us now change the order of the `rspde_model`

object to
be 1:

```
rational.order(rspde_model) <- 1
rational.order(Abar) <- 1
rational.order(mesh.index) <- 1
```

Let us fit this new model:

```
.1 <- y ~ -1 + Intercept + f(seaDist, model = "rw1") +
f.sf(field, model = rspde_model)
.1 <- inla.stack(
stk.datdata = list(y = Y), A = list(Abar, 1), tag = "est",
effects = list(
c(
mesh.index,list(Intercept = 1)
),list(
long = inla.group(coords[, 1]),
lat = inla.group(coords[, 2]),
seaDist = inla.group(seaDist)
)
)
)
<- inla(f.s.1,
rspde_fit_order_1 family = "Gamma", data = inla.stack.data(stk.dat.1),
verbose = FALSE,
control.inla = list(int.strategy = "eb"),
control.predictor = list(A = inla.stack.A(stk.dat.1), compute = TRUE)
)
```

Here is the summary:

`summary(rspde_fit_order_1)`

```
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 3.17, Running = 5.68, Post = 0.0521, Total = 8.91
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 0.972 0.026 0.92 0.972 1.023 0.972 0
##
## Random effects:
## Name Model
## seaDist RW1 model
## field CGeneric
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision parameter for the Gamma observations 13.220 0.901 11.491
## Precision for seaDist 9712.216 7301.457 2427.830
## Theta1 for field -1.478 0.147 -1.777
## Theta2 for field -0.054 0.300 -0.631
## Theta3 for field -0.301 0.890 -1.999
## 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 13.205 15.040 13.205
## Precision for seaDist 7693.402 29111.618 5131.661
## Theta1 for field -1.475 -1.196 -1.463
## Theta2 for field -0.059 0.551 -0.078
## Theta3 for field -0.323 1.524 -0.416
##
## Marginal log-Likelihood: -1263.10
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
```

We can fix the smoothness, say \(\nu\), of the model by providing a non-NULL
positive value for `nu`

.

When the smoothness, \(\nu\), is fixed, we can have two possibilities:

\(\alpha = \nu + d/2\) is integer;

\(\alpha = \nu + d/2\) is not integer.

The first case, i.e., when \(\alpha\) is integer, has less computational cost. Furthermore, the \(A\) matrix is different than the \(A\) matrix for the non-integer \(\alpha\).

The \(A\) matrix is the same for all
values of \(\nu\) such that \(\alpha\) is integer. So, the \(A\) matrix for these cases only need to be
computed once. The same holds for the `index`

obtained from
the `rspde.make.index()`

function.

In the second case the \(A\) matrix
only depends on the order of the rational approximation and not on \(\nu\). Therefore, if the matrix \(A\) has already been computed for some
`rspde.order`

, then the \(A\) matrix will be same for all the values
of \(\nu\) such that \(\alpha\) is non-integer for that
`rspde.order`

. The same holds for the `index`

obtained from the `rspde.make.index()`

function.

If \(\nu\) is fixed, we have that
the parameters returned by `R-INLA`

are \[\kappa = \exp(\theta_1)\quad\hbox{and}\quad\tau =
\exp(\theta_2).\] We will now provide illustrations for both
scenarios. It is also noteworthy that just as for the case in which we
estimate \(\nu\), we can also change
the order of the rational approximation by changing the value of
`rspde.order`

. For both illustrations with fixed \(\nu\) below, we will only consider the
order of the rational approximation of 2, that is, the default
order.

Recall that: \[\nu =
\nu_{UB}\Big(\frac{\exp(\theta_3)}{1+\exp(\theta_3)}\Big).\]
Thus, to illustrate, let us consider a fixed \(\nu\) given by the mean of \(\nu\) obtained from the first model we
considered in this vignette, namely, the fit given by
`rspde_fit`

, which is approximately \(\nu = 1.21\).

Notice that for this \(\nu\), the value of \(\alpha\) is non-integer, so we can use the \(A\) matrix and the index of the first fitted model, which is also of order 2.

Therefore, all we have to do is build a new model in which we set
`nu`

to `1.21`

:

```
<- rspde.matern(mesh = prmesh, rspde.order = 2,
rspde_model_fix nu = 1.21
)
```

Let us now fit the model:

```
<- y ~ -1 + Intercept + f(seaDist, model = "rw1") +
f.s.fix f(field, model = rspde_model_fix)
<- inla(f.s.fix,
rspde_fix family = "Gamma", data = inla.stack.data(stk.dat),
verbose = FALSE,
control.inla = list(int.strategy = "eb"),
control.predictor = list(A = inla.stack.A(stk.dat), compute = TRUE)
)
```

Here we have the summary:

`summary(rspde_fix)`

```
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 3.17, Running = 4.12, Post = 0.052, Total = 7.35
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 0.648 0.018 0.614 0.648 0.682 0.648 0
##
## Random effects:
## Name Model
## seaDist RW1 model
## field CGeneric
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision parameter for the Gamma observations 13.186 0.898 11.466
## Precision for seaDist 9978.781 7897.570 2468.399
## Theta1 for field -1.480 0.140 -1.756
## Theta2 for field 0.004 0.296 -0.545
## 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 13.172 15.003 13.168
## Precision for seaDist 7737.047 30970.995 5064.285
## Theta1 for field -1.480 -1.204 -1.480
## Theta2 for field -0.008 0.618 -0.056
##
## Marginal log-Likelihood: -1263.73
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
```

Now, the summary in the original scale:

```
<- rspde.result(rspde_fix, "field", rspde_model_fix)
result_fix summary(result_fix)
```

```
## mean sd 0.025quant 0.5quant 0.975quant mode
## std.dev 0.229865 0.0321689 0.173121 0.227668 0.29926 0.223471
## range 1.048880 0.3242590 0.583023 0.989949 1.84336 0.883016
```

Since we are in dimension \(d=2\), and \(\nu>0\), the smallest value of \(\nu\) that makes \(\alpha = \nu + 1\) an integer is \(\nu=1\). This value is also close to the estimated mean of the first model we fitted and enclosed by the posterior marginal density of \(\nu\) for the first fit. Therefore, let us fit the model with \(\nu=1\).

To this end we need to compute a new \(A\) matrix:

```
<- rspde.make.A(
Abar.int mesh = prmesh, loc = coords,
nu = 1
)
```

a new index:

```
<- rspde.make.index(
mesh.index.int name = "field", mesh = prmesh,
nu = 1
)
```

create a new model (remember to set `nu=1`

):

```
<- rspde.matern(mesh = prmesh,
rspde_model_fix_int1 nu = 1)
```

The remaining is standard:

```
<- inla.stack(
stk.dat.int data = list(y = Y), A = list(Abar.int, 1), tag = "est",
effects = list(
c(
mesh.index.int,list(Intercept = 1)
),list(
long = inla.group(coords[, 1]),
lat = inla.group(coords[, 2]),
seaDist = inla.group(seaDist)
)
)
)
.1 <- y ~ -1 + Intercept + f(seaDist, model = "rw1") +
f.s.fix.intf(field, model = rspde_model_fix_int1)
<- inla(f.s.fix.int.1,
rspde_fix_int_1 family = "Gamma",
data = inla.stack.data(stk.dat.int), verbose = FALSE,
control.inla = list(int.strategy = "eb"),
control.predictor = list(
A = inla.stack.A(stk.dat.int),
compute = TRUE
) )
```

Let us check the summary:

`summary(rspde_fix_int_1)`

```
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 3.22, Running = 1.14, Post = 0.0407, Total = 4.4
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 1.944 0.055 1.837 1.944 2.052 1.944 0
##
## Random effects:
## Name Model
## seaDist RW1 model
## field CGeneric
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision parameter for the Gamma observations 13.186 0.896 11.487
## Precision for seaDist 10949.795 9956.683 2502.565
## Theta1 for field -1.457 0.142 -1.735
## Theta2 for field 0.056 0.326 -0.543
## 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 13.16 15.017 13.135
## Precision for seaDist 8007.49 37218.662 4942.591
## Theta1 for field -1.46 -1.176 -1.459
## Theta2 for field 0.04 0.738 -0.021
##
## Marginal log-Likelihood: -1262.45
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
```

and check the summary in the user’s scale:

```
<- rspde.result(rspde_fix_int_1, "field", rspde_model_fix_int1)
rspde_result_int summary(rspde_result_int)
```

```
## mean sd 0.025quant 0.5quant 0.975quant mode
## std.dev 0.235375 0.0333692 0.176869 0.232936 0.307672 0.228078
## range 1.115750 0.3855130 0.584329 1.038740 2.077980 0.901474
```

We begin by recalling that the fitted `rSPDE`

model in `R-INLA`

contains the
parameters \(\textrm{Theta1}\), \(\textrm{Theta2}\) and \(\textrm{Theta3}\). Let, again, \(\theta_1 = \textrm{Theta1}\), \(\theta_2=\textrm{Theta2}\) and \(\theta_3=\textrm{Theta3}\). In terms of the
SPDE \[(\kappa^2 I - \Delta)^{\alpha/2}(\tau
u) = \mathcal{W},\] where \(\alpha =
\nu + d/2\).

We also have the range parameter \(\rho = \frac{\sqrt{8\nu}}{\kappa}\) and the standard deviation \(\sigma = \sqrt{\frac{\Gamma(\nu)}{\tau^2 \kappa^{2\nu}(4\pi)^{d/2}\Gamma(\nu + d/2)}}\).

We begin by dealing with \(\tau\) and \(\kappa\).

We have that \[\tau = \exp(\theta_1),\quad
\kappa = \exp(\theta_2).\] The `rspde.matern()`

function assumes a lognormal prior distribution for \(\tau\) and \(\kappa\). This prior distribution is
obtained by assuming that \(\theta_1\)
and \(\theta_2\) follow normal
distributions. By default we assume \(\theta_1\) and \(\theta_2\) to be independent and to follow
normal distributions \(\theta_1\sim
N(\log(\tau_0), 10)\) and \(\theta_2\sim N(\log(\kappa_0), 10)\).

\(\kappa_0\) is suitably defined in terms of the mesh and \(\tau_0\) is defined in terms of \(\kappa_0\) and on the prior of the smoothness parameter.

If one wants to define the prior \[\theta_1 \sim N(\text{mean_theta_1},
\text{sd_theta_1}),\] one can simply set the argument
`prior.tau = list(meanlog=mean_theta_1, sdlog=sd_theta_1)`

.
Analogously, to define the prior \[\theta_2
\sim N(\text{mean_theta_2}, \text{sd_theta_2}),\] one can set the
argument
`prior.kappa = list(meanlog=mean_theta_2, sdlog=sd_theta_2)`

.

It is important to mention that, by default, the initial values of \(\tau\) and \(\kappa\) are \(\exp(\text{mean_theta_1})\) and \(\exp(\text{mean_theta_2})\), respectively. So, if the user does not change these parameters, and also does not change the initial values, the initial values of \(\tau\) and \(\kappa\) will be, respectively, \(\tau_0\) and \(\kappa_0\).

If one sets `prior.tau = list(meanlog=mean_theta_1)`

, the
prior for \(\theta_1\) will be \[\theta_1 \sim N(\text{mean_theta_1}, 1),\]
whereas, if one sets, `prior.tau = list(sdlog=sd_theta_1)`

,
the prior will be \[\theta_1 \sim
N(\log(\tau_0), \text{sd_theta_1}).\] Analogously, if one sets
`prior.kappa = list(meanlog=mean_theta_2)`

, the prior for
\(\theta_2\) will be \[\theta_2 \sim N(\text{mean_theta_2}, 1),\]
whereas, if one sets, `prior.kappa = list(sdlog=sd_theta_2)`

,
the prior will be \[\theta_2 \sim
N(\log(\kappa_0), \text{sd_theta_2}).\]

Let us now consider the priors for the range, \(\rho\), and std. deviation, \(\sigma\). This parameterization is used
with the argument `parameterization = "matern"`

, which is the
default.

In this case, we have that \[\sigma =
\exp(\theta_1),\quad \rho = \exp(\theta_2).\] We have two options
for prior. By default, which is the option
`prior.theta.param = "theta"`

, the
`rspde.matern()`

function assumes a lognormal prior
distribution for \(\sigma\) and \(\rho\). This prior distribution is obtained
by assuming that \(\theta_1\) and \(\theta_2\) follow normal distributions. By
default we assume \(\theta_1\) and
\(\theta_2\) to be independent and to
follow normal distributions \(\theta_1\sim
N(\log(\sigma_0), 10)\) and \(\theta_2\sim N(\log(\rho_0), 10)\).

\(\rho_0\) is suitably defined in terms of the mesh and \(\sigma_0\) is defined in terms of \(\rho_0\) and on the prior of the smoothness parameter.

Similarly to the priors of \(\tau\)
and \(\kappa\), if one wants to define
the prior \[\theta_1 \sim
N(\text{mean_theta_1}, \text{sd_theta_1}),\] one can simply set
the argument
`prior.tau = list(meanlog=mean_theta_1, sdlog=sd_theta_1)`

.
Analogously, to define the prior \[\theta_2
\sim N(\text{mean_theta_2}, \text{sd_theta_2}),\] one can set the
argument
`prior.kappa = list(meanlog=mean_theta_2, sdlog=sd_theta_2)`

.

Another option is to set `prior.theta.param = "spde"`

. In
this case, a change of variables is performed. So, we assume a lognormal
prior on \(\tau\) and \(\kappa\). Then, by the relations \(\rho = \frac{\sqrt{8\nu}}{\kappa}\) and
\(\sigma = \sqrt{\frac{\Gamma(\nu)}{\tau^2
\kappa^{2\nu}(4\pi)^{d/2}\Gamma(\nu + d/2)}}\), we obtain a prior
for \(\rho\) and \(\sigma\).

Finally, let us consider the smoothness parameter \(\nu\).

By default, we assume that \(\nu\) follows a truncated lognormal distribution. The truncated lognormal distribution is defined in the following sense. We assume that \(\log(\nu)\) has prior distribution given by a truncated normal distribution with support \((-\infty,\log(\nu_{UB}))\), where \(\nu_{UB}\) is the upper bound for \(\nu\), with location parameter \(\mu_0 =\log(\nu_0)= \log\Big(\min\{1,\nu_{UB}/2\}\Big)\) and scale parameter \(\sigma_0 = 1\). More precisely, let \(\Phi(\cdot; \mu,\sigma)\) stand for the cumulative distribution function (CDF) of a normal distribution with mean \(\mu\) and standard deviation \(\sigma\). Then, \(\log(\nu)\) has cumulative distribution function given by \[F_{\log(\nu)}(x) = \frac{\Phi(x;\mu_0,\sigma_0)}{\Phi(\nu_{UB})},\quad x\leq \nu_{UB},\] and \(F_{\log(\nu)}(x) = 1\) if \(x>\nu_{UB}\). We will call \(\mu_0\) and \(\sigma_0\) the log-location and log-scale parameters of \(\nu\), respectively, and we say that \(\log(\nu)\) follows a truncated normal distribution with location parameter \(\mu_0\) and scale parameter \(\sigma_0\).

We can have another possibility of prior distribution for \(\nu\), namely, a beta distribution on the
interval \((0,\nu_{UB})\), where \(\nu_{UB}\) is the upper bound for \(\nu\), with mean \(\nu_0=\min\{1, \nu_{UB}/2\}\) and variance
\(\frac{\nu_0(\nu_{UB}-\nu_0)}{1+\phi_0}\),
and we call \(\phi_0\) the precision
parameter, whose default value is \[\phi_0 =
\max\Big\{\frac{\nu_{UB}}{\nu_0}, \frac{\nu_{UB}}{\nu_{UB}-\nu_0}\Big\}
+ \phi_{inc}.\] The parameter \(\phi_{inc}\) is an increment to ensure that
the prior beta density has boundary values equal to zero (where the
boundary values are defined either by continuity or by limits). The
default value of \(\phi_{inc}\) is 1.
The value of \(\phi_{inc}\) can be
changed by changing the argument `nu.prec.inc`

in the
`rspde.matern()`

function. The higher the value of \(\phi_{inc}\) (that is, the value of
`nu.prec.inc`

) the more informative the prior distribution
becomes.

Let us denote a beta distribution with support on \((0,\nu_{UB})\), mean \(\mu\) and precision parameter \(\phi\) by \(\mathcal{B}_{\nu_{UB}}(\mu,\phi)\).

If we want \(\nu\) to have a prior
\[\nu \sim
\mathcal{B}_{\nu_{UB}}(\text{nu_1},\text{prec_1}),\] one simply
needs to set `prior.nu = list(mean=nu_1, prec=prec_1)`

. If
one sets `prior.nu = list(mean=nu_1)`

, then \(\nu\) will have prior \[\nu \sim
\mathcal{B}_{\nu_{UB}}(\text{nu_1},\phi_1),\] where \[\phi_1 = \max\Big\{\frac{\nu_{UB}}{\text{nu_1}},
\frac{\nu_{UB}}{\nu_{UB}-\text{nu_1}}\Big\} +
\text{nu.prec.inc}.\]

Of one sets `prior.nu = list(prec=prec_1)`

, then \(\nu\) will have prior \[\nu\sim \mathcal{B}_{\nu_{UB}}(\nu_0,
\text{prec_1}).\] It is also noteworthy that we have that, in
terms of `R-INLA`

’s
parameters,

\[\nu = \nu_{UB}\Big(\frac{\exp(\theta_3)}{1+\exp(\theta_3)}\Big).\]

It is important to mention that, by default, if a beta prior distribution is chosen for the smoothness parameter \(\nu\), then the initial value of \(\nu\) is the mean of the prior beta distribution. So, if the user does not change this parameter, and also does not change the initial value, the initial value of \(\nu\) will be \(\min\{1,\nu_{UB}/2\}\).

We also assume that, in terms of `R-INLA`

’s parameters,
\[\nu =
\nu_{UB}\Big(\frac{\exp(\theta_3)}{1+\exp(\theta_3)}\Big).\]

To change the prior distribution of \(\nu\) to the truncated lognormal
distribution, we need to set the argument
`prior.nu.dist="lognormal"`

.

To change these parameters in the prior distribution to, say,
`log_nu_1`

and `log_sigma_1`

, one can simply set
`prior.nu = list(loglocation=log_nu_1, logscale=sigma_1)`

.

If one sets `prior.nu = list(loglocation=log_nu_1)`

, the
prior for \(\theta_3\) will be a
truncated normal normal distribution with location parameter
`log_nu_1`

and scale parameter `1`

. Analogously,
if one sets, `prior.nu = list(logscale=sigma_1)`

, the prior
for \(\theta_3\) will be a truncated
normal distribution with location parameter \(\log(\nu_0)=
\log\Big(\min\{1,\nu_{UB}/2\}\Big)\) and scale parameter
`sigma_1`

.

It is important to mention that, by default, if a truncated lognormal prior distribution is chosen for the smoothness parameter \(\nu\), then the initial value of \(\nu\) is the exponential of the log-location parameter of \(\nu\). So, if the user does not change this parameter, and also does not change the initial value, the initial value of \(\nu\) will be \(\min\{1,\nu_{UB}/2\}\).

Let us consider an example with the same dataset used in the first model of this vignette where we change the prior distribution of \(\nu\) from beta to lognormal.

`<- rspde.matern(mesh = prmesh, prior.nu.dist = "beta") rspde_model_beta `

Since we did not change `rspde.order`

and are not fixing
\(\nu\), we can use the same \(A\) matrix and same index from the first
example.

Therefore, all we have to do is update the formula and fit the model:

```
<- y ~ -1 + Intercept + f(seaDist, model = "rw1") +
f.s.beta f(field, model = rspde_model_beta)
<- inla(f.s.beta,
rspde_fit_beta family = "Gamma", data = inla.stack.data(stk.dat),
verbose = FALSE,
control.inla = list(int.strategy = "eb"),
control.predictor = list(A = inla.stack.A(stk.dat), compute = TRUE)
)
```

We have the summary:

`summary(rspde_fit_beta)`

```
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 3.24, Running = 17.5, Post = 0.0588, Total = 20.8
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 0.639 0.358 -0.063 0.639 1.342 0.639 0
##
## Random effects:
## Name Model
## seaDist RW1 model
## field CGeneric
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision parameter for the Gamma observations 13.440 0.917 11.73
## Precision for seaDist 22453.925 9675.038 7660.86
## Theta1 for field -0.863 0.059 -1.01
## Theta2 for field 2.422 0.309 1.66
## Theta3 for field -2.621 0.267 -3.28
## 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 13.405 15.35 13.329
## Precision for seaDist 20920.358 45533.80 17895.016
## Theta1 for field -0.879 -0.74 -0.845
## Theta2 for field 2.316 3.08 2.513
## Theta3 for field -2.697 -2.06 -2.540
##
## Marginal log-Likelihood: -1258.70
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
```

Also, we can have the summary in the user’s scale:

```
<- rspde.result(rspde_fit_beta, "field", rspde_model_beta)
result_fit_beta summary(result_fit_beta)
```

```
## mean sd 0.025quant 0.5quant 0.975quant mode
## std.dev 0.422763 0.0243023 0.371497 0.424389 0.466299 0.428050
## range 11.796000 3.4437600 5.754750 11.605900 19.094200 11.326400
## nu 0.209193 0.0493008 0.117411 0.208405 0.308472 0.209176
```

and the plot of the posterior marginal densities

```
<- gg_df(result_fit_beta)
posterior_df_fit_beta
ggplot(posterior_df_fit_beta) + geom_line(aes(x = x, y = y)) +
facet_wrap(~parameter, scales = "free") + labs(y = "Density")
```

The starting values to be used by `R-INLA`

’s optimization
algorithm can be changed by setting the arguments
`start.ltau`

, `start.lkappa`

and
`start.nu`

.

`start.ltau`

will be the initial value for \(\log(\tau)\), that is, the logarithm of \(\tau\).`start.lkappa`

will be the inital value for \(\log(\kappa)\), that is, the logarithm of \(\kappa\).`start.nu`

will be the initial value for \(\nu\). Notice that here the initial value is*not*on the log scale.

One can change the initial value of one or more parameters.

For instance, let us consider the example with precipitation data,
`rspde.order=3`

, but change the initial values to the ones
close to the fitted value when considering the default
`rspde.order`

(which is 2):

```
<- rspde.matern(mesh = prmesh, rspde.order = 3,
rspde_model_order_3_start nu.upper.bound = 2,
start.lkappa = result_fit$summary.log.kappa$mean,
start.ltau = result_fit$summary.log.tau$mean,
start.nu = min(result_fit$summary.nu$mean, 2 - 1e-5)
)
```

Since we already computed the \(A\)
matrix and the index for `rspde.order=3`

, all we have to do
is to update the formula and fit:

```
3.start <- y ~ -1 + Intercept + f(seaDist, model = "rw1") +
f.s.f(field, model = rspde_model_order_3_start)
<- inla(f.s.3.start,
rspde_fit_order_3_start family = "Gamma",
data = inla.stack.data(stk.dat.3),
verbose = FALSE,
control.inla = list(int.strategy = "eb"),
control.predictor = list(
A = inla.stack.A(stk.dat.3),
compute = TRUE
) )
```

We have the summary:

`summary(rspde_fit_order_3_start)`

```
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 3.18, Running = 6.69, Post = 0.0633, Total = 9.94
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 0.486 0.014 0.458 0.486 0.514 0.486 0
##
## Random effects:
## Name Model
## seaDist RW1 model
## field CGeneric
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision parameter for the Gamma observations 13.189 0.896 11.494
## Precision for seaDist 11204.965 9615.651 2553.217
## Theta1 for field -1.430 0.162 -1.742
## Theta2 for field 0.103 0.393 -0.626
## Theta3 for field -0.486 0.854 -2.208
## 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 13.166 15.023 13.132
## Precision for seaDist 8425.761 36773.944 5289.497
## Theta1 for field -1.432 -1.102 -1.443
## Theta2 for field 0.087 0.921 0.022
## Theta3 for field -0.473 1.166 -0.423
##
## Marginal log-Likelihood: -1263.16
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
```

We have three rational approximations available. The BRASIL algorithm Hofreither (2021), and two “versions” of the Clenshaw-Lord Chebyshev-Pade algorithm, one with lower bound zero and another with the lower bound given in Xiong, Simas, and Bolin (2022).

The type of rational approximation can be chosen by setting the
`type.rational.approx`

argument in the
`rspde.matern`

function. The BRASIL algorithm corresponds to
the choice `brasil`

, the Clenshaw-Lord Chebyshev pade with
zero lower bound and non-zero lower bounds are given, respectively, by
the choices `chebfun`

and `chebfunLB`

.

Let us fit a model assigning a `brasil`

rational
approximation. We will consider a model with the order of the rational
approximation being 1:

```
<- rspde.matern(prmesh,
rspde_model_brasil type.rational.approx = "brasil")
<- y ~ -1 + Intercept + f(seaDist, model = "rw1") +
f.s.brasil f(field, model = rspde_model_brasil)
<- inla(f.s.brasil,
rspde_fit_order_1_brasil family = "Gamma", data = inla.stack.data(stk.dat),
verbose = FALSE,
control.inla = list(int.strategy = "eb"),
control.predictor = list(A = inla.stack.A(stk.dat), compute = TRUE)
)
```

Let us get the summary:

`summary(rspde_fit_order_1_brasil)`

```
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 3.22, Running = 12.9, Post = 0.0559, Total = 16.2
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 0.648 0.019 0.611 0.648 0.685 0.648 0
##
## Random effects:
## Name Model
## seaDist RW1 model
## field CGeneric
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision parameter for the Gamma observations 13.165 0.898 11.446
## Precision for seaDist 10056.140 8598.481 2369.512
## Theta1 for field -1.461 0.178 -1.809
## Theta2 for field 0.048 0.366 -0.656
## Theta3 for field -0.697 0.918 -2.511
## 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 13.150 14.981 13.147
## Precision for seaDist 7551.322 32842.899 4779.909
## Theta1 for field -1.461 -1.110 -1.463
## Theta2 for field 0.042 0.785 0.018
## Theta3 for field -0.694 1.103 -0.684
##
## Marginal log-Likelihood: -1263.90
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
```

Finally, similarly to the order of the rational approximation, one
can check the order with the `rational.type()`

function, and
assign a new type with the `rational.type<-()`

function.

`rational.type(rspde_model)`

`## [1] "chebfun"`

`rational.type(rspde_model_brasil)`

`## [1] "brasil"`

Let us change the type of the rational approximation on the model with rational approximation of order 3:

```
rational.type(rspde_model_order_3) <- "brasil"
.3 <- y ~ -1 + Intercept + f(seaDist, model = "rw1") +
f.sf(field, model = rspde_model_order_3)
<- inla(f.s.3,
rspde_fit_order_3_brasil family = "Gamma", data = inla.stack.data(stk.dat.3),
verbose = FALSE,
control.inla = list(int.strategy = "eb"),
control.predictor = list(A = inla.stack.A(stk.dat.3), compute = TRUE)
)
```

Let us get the summary:

`summary(rspde_fit_order_3_brasil)`

```
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 3.16, Running = 9.07, Post = 0.089, Total = 12.3
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 0.486 0.014 0.458 0.486 0.514 0.486 0
##
## Random effects:
## Name Model
## seaDist RW1 model
## field CGeneric
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision parameter for the Gamma observations 13.234 0.902 11.509
## Precision for seaDist 10848.797 8207.669 2858.223
## Theta1 for field -1.413 0.152 -1.707
## Theta2 for field 0.100 0.400 -0.647
## Theta3 for field -0.594 0.788 -2.196
## 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 13.219 15.060 13.213
## Precision for seaDist 8555.274 32765.188 5724.690
## Theta1 for field -1.415 -1.106 -1.425
## Theta2 for field 0.085 0.934 0.022
## Theta3 for field -0.581 0.934 -0.530
##
## Marginal log-Likelihood: -1263.02
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
```

Hofreither, Clemens. 2021. “An Algorithm for Best Rational
Approximation Based on Barycentric Rational Interpolation.”
*Numerical Algorithms* 88 (1): 365–88.

Lindgren, Finn, Håvard Rue, and Johan 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–98.

Xiong, Zhen, Alexandre B. Simas, and David Bolin. 2022.
“Covariance-Based Rational Approximations of Fractional SPDEs for
Computationally Efficient Bayesian Inference.”
*arXiv:2209.04670*.