In this vignette we provide a brief introduction to the `rSPDE`

package. We begin by using the package to generate a simple data set. Then, we will analyze this data set using the `rSPDE`

package.

We have two main “families” of functions inside the `rSPDE`

package:

The

`R-INLA`

implementation of the rational SPDE approach of (Bolin and Kirchner 2020)(https://www.tandfonline.com/doi/full/10.1080/10618600.2019.1665537);The rational SPDE approach.

We will analyze the same data set using both approaches. We also have specific vignettes for each of the aforementioned family. More precisely, we have the following additional vignettes:

We also have a separate group of functions for performing operator-based rational approximations. These will be useful when performing rational approximations for fractional SPDE models with non-Gaussian noise. An example in which such approximation is suitable is when one has the so-called type-G Lévy noises.

We refer the reader to (Wallin and Bolin 2015)(https://onlinelibrary.wiley.com/doi/full/10.1111/sjos.12141), (Bolin 2013)(https://onlinelibrary.wiley.com/doi/abs/10.1111/sjos.12046) and (Asar et al. 2020)(https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/rssc.12405) for examples of models driven by type-G Lévy noises. We also refer the reader to the `ngme`

package where one can fit such models.

We explore the functions for performing the operator-based rational approximation on the vignette:

We begin by generating a toy data set.

For this illustration, we will simulate a data set on a two-dimensional spatial domain. To this end, we need to construct a mesh over the domain of interest and then compute the matrices needed to define the operator. We will use the `R-INLA`

package to create the mesh and obtain the matrices of interest.

We will begin by defining a mesh over \([0,1]\times [0, 1]\):

```
library(INLA)
m = 500
loc_2d_mesh = matrix(runif(m*2),m,2)
mesh_2d = inla.mesh.2d(
loc=loc_2d_mesh,
cutoff=0.05,
max.edge=c(0.1,0.5) )
plot(mesh_2d, main = "")
```

We will now use the `matern.operators()`

function to construct a rational SPDE approximation of order \(m=2\) for a Gaussian random field with a Matérn covariance function on \([0,1]\times [0, 1]\). We choose \(\nu=0.5\) which corresponds to exponential covariance. We also set \(\sigma=1\) and the range as \(0.2\).

```
library(rSPDE)
sigma <- 1
range <- 0.2
nu <- 0.5
kappa <- sqrt(8*nu)/range
op <- matern.operators(mesh=mesh_2d,nu=nu,
kappa=kappa,sigma=sigma,m=2)
```

We are now in a position to simulate the latent field \(u\):

Let us then consider a simple Gaussian linear model of a latent spatial field \(u(\mathbf{s})\), observed at \(m\) locations, \(\{\mathbf{s}_1 , \ldots , \mathbf{s}_m \}\). For each \(i = 1,\ldots,m,\) we have \[ \begin{align} y_i &= u(\mathbf{s}_i)+\varepsilon_i\\ \end{align}, \] where \(\varepsilon_1,\ldots,\varepsilon_{m}\) are iid normally distributed with mean 0 and standard deviation 0.1.

We now obtain a realization of \(y\). We will use the `R-INLA`

function `inla.spde.make.A()`

to construct the observation matrix.

```
A <- inla.spde.make.A(
mesh=mesh_2d,
loc=loc_2d_mesh)
sigma.e <- 0.1
y = A %*% u + rnorm(m) * sigma.e
```

The observed data can be seen in the following image.

The simulated random field is shown in the following figure.

`R-INLA`

implementation of the rational SPDE approachWe will now fit the model of the toy data set using our `R-INLA`

implementation of the rational SPDE approach. Further details on this implementation can be found in R-INLA implementation of the rational SPDE approach.

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

object.

```
Abar <- rspde.make.A(mesh = mesh_2d, loc = loc_2d_mesh)
mesh.index <- rspde.make.index(name = "field", mesh = mesh_2d)
st.dat=inla.stack(
data=list(y=as.vector(y)),
A=Abar,
effects=mesh.index)
```

We now create the model object. We need to set an upper bound for the smoothness parameter \(\nu\). When we increase the upper bound for \(\nu\) we also increase the computational cost. For this example we set `nu_upper_bound=1`

. See the R-INLA implementation of the rational SPDE approach for further details.

Finally, we create the formula and fit:

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

We can get the summary:

```
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 = 4.32, Running = 46.9, Post = 0.0831, Total = 51.3
#> Random effects:
#> Name Model
#> field RGeneric2
#>
#> Model hyperparameters:
#> mean sd 0.025quant 0.5quant
#> Precision for the Gaussian observations 75.623 3.400 67.831 76.173
#> Theta1 for field -1.972 0.068 -2.086 -1.979
#> Theta2 for field 2.454 0.021 2.411 2.454
#> Theta3 for field 0.709 0.020 0.669 0.709
#> 0.975quant mode
#> Precision for the Gaussian observations 80.50 78.856
#> Theta1 for field -1.82 -2.007
#> Theta2 for field 2.50 2.454
#> Theta3 for field 0.75 0.709
#>
#> Marginal log-Likelihood: 28.81
#> 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:

```
result_fit <- rspde.result(rspde_fit, "field", rspde_model)
summary(result_fit)
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> tau 0.139349 0.00961679 0.124223 0.138200 0.161752 0.134029
#> kappa 11.626600 0.24451800 11.166800 11.636100 12.130400 11.629900
#> nu 0.669521 0.00451107 0.661483 0.670229 0.679032 0.670204
tau <- op$tau
result_df <- data.frame(parameter = c("tau","kappa","nu"),
true = c(tau, kappa, nu), mean = c(result_fit$summary.tau$mean,
result_fit$summary.kappa$mean,
result_fit$summary.nu$mean),
mode = c(result_fit$summary.tau$mode,
result_fit$summary.kappa$mode,
result_fit$summary.nu$mode))
print(result_df)
#> parameter true mean mode
#> 1 tau 0.2236068 0.1393486 0.1340293
#> 2 kappa 10.0000000 11.6266083 11.6298835
#> 3 nu 0.5000000 0.6695206 0.6702036
```

`R-INLA`

implementation of the rational SPDE approachLet us now obtain predictions (i.e., do kriging) of the latent field 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()`

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:

This lattice contains 100 × 100 locations (the default).

Let us plot the locations that we will do prediction:

Let us now calculate the predictions jointly with the estimation. To this end, first, we begin by linking the prediction coordinates to the mesh nodes through an \(A\) matrix

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.

```
ef.prd = list(c(mesh.index))
st.prd <- inla.stack(data = list(y = NA),
A = list(A.prd), tag = "prd",
effects = ef.prd)
st.all <- inla.stack(st.dat, st.prd)
```

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:

```
rspde_fitprd <- inla(f, family = "Gaussian",
data = inla.stack.data(st.all),
control.predictor = list(A = inla.stack.A(st.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:

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

Finally, we plot the results. First the mean:

Then, the std. deviations:

`rSPDE`

We will now fit the model of the toy data set without using `R-INLA`

. To this end we will use the rational approximation functions from `rSPDE`

package. Further details can be found in the vignette Rational approximation with the rSPDE package.

We will now use the function `rSPDE.matern.loglike()`

to define the likelihood. This function is object-based, in the sense that it obtains several of the quantities it needs from the `rSPDE`

model object.

Notice that we already created a `rSPDE`

model object to simulate the data. We will, then, use the same model object. Recall that the `rSPDE`

model object we created is `op`

. We also already have the \(A\) matrix connecting the observation locations to the mesh, and we simply called it `A`

.

To simplify parameter estimation, we create an objective function to minimize which is the negative log-likelihood, parametrized using the logarithm of each parameter to avoid constrained optimization.

```
mlik <- function(theta, Y, A, op) {
sigma = exp(theta[1])
kappa = exp(theta[2])
nu = exp(theta[3])
return(-rSPDE::rSPDE.matern.loglike(object = op, Y=Y,
A = A, user_kappa=kappa, user_sigma = sigma,
user_nu=nu, sigma.e = exp(theta[4])))
}
```

We can now estimate the parameter using `optimParallel()`

(one can also use `optim()`

):

```
library(optimParallel)
#Preparing the parallel
#Checking if we have a limit to the number of cores
chk <- Sys.getenv("_R_CHECK_LIMIT_CORES_", "")
if (nzchar(chk) && chk == "TRUE") {
n_cores <- 2L
} else {
n_cores <- parallel::detectCores() - 1
}
cl <- makeCluster(n_cores)
setDefaultCluster(cl=cl)
#Fitting the model
theta0 = c(get.inital.values.rSPDE(mesh=mesh_2d),
log(0.1*sqrt(var(as.vector(y)))))
start_time <- Sys.time()
pars <- optimParallel(theta0, mlik, Y = y, A = A, op=op)
end_time <- Sys.time()
total_time <- end_time - start_time
results <- data.frame(sigma = c(sigma, exp(pars$par[1])),
kappa = c(kappa, exp(pars$par[2])),
nu = c(nu, exp(pars$par[3])),
sigma.e = c(sigma.e, exp(pars$par[4])),
row.names = c("True", "Estimate"))
print(results)
#> sigma kappa nu sigma.e
#> True 1 10.000000 0.5000000 0.1000000
#> Estimate 1 8.098931 0.5277972 0.1134149
#Total time
print(total_time)
#> Time difference of 7.262474 secs
```

`rSPDE`

We will now do kringing on the same dense grid we did for the `R-INLA`

-based rational SPDE approach, but now using the `rSPDE`

functions. To this end we will use the `predict`

method on the `rSPDE`

model object.

Observe that we need an \(A\) matrix connecting the mesh to the prediction locations.

Let us now create the \(A\) matrix for the same prediction locations we used for the previous case (using the `R-INLA`

implementation):

We will now use the `predict()`

method on the `rSPDE`

model object with the argument `compute.variances`

set to `TRUE`

so that we can plot the standard deviations. Let us also update the values of the `rSPDE`

model object to the fitted ones, and also save the estimated value of `sigma.e`

.

```
sigma.e.est <- exp(pars$par[4])
op.prd <- update(op, user_sigma = exp(pars$par[1]),
user_kappa = exp(pars$par[2]),
user_nu = exp(pars$par[3]))
pred.rspde <- predict(op.prd, A = A, Aprd = A.prd2, Y = y,
sigma.e = sigma.e.est,
compute.variances = TRUE)
```

Finally, we plot the results. First the mean:

Then, the std. deviations:

Asar, Özgür, David Bolin, Peter Diggle, and Jonas Wallin. 2020. “Linear Mixed Effects Models for Non‐Gaussian Repeated Measurement Data.” *Journal of the Royal Statistical Society. Series C. Applied Statistics* 69 (5): 1015–65.

Bolin, David. 2013. “Spatial Matérn Fields Driven by Non-Gaussian Noise.” *Scandinavian Journal of Statistics* 41 (3): 557–79.

Bolin, David, and Kristin Kirchner. 2020. “The Rational SPDE Approach for Gaussian Random Fields with General Smoothness.” *Journal of Computational and Graphical Statistics* 29 (2): 274–85.

Wallin, Jonas, and David Bolin. 2015. “Geostatistical Modelling Using Non-Gaussian Matérn Fields.” *Scandinavian Journal of Statistics* 42 (3): 872–90.