An introduction to the rSPDE package

David Bolin and Alexandre B. Simas

2023-06-25

Introduction

In this vignette we provide a brief introduction to the rSPDE package. The main approach for constructing the rational approximations is the covariance-based rational SPDE approach of Bolin, Simas, and Xiong (2023). The package contains three main “families” of functions that implement the approach:

To illustrate these different functions, we begin by using the package to generate a simple data set, which then will be analyzed using the different approaches. Further details on each family of functions is given in the following additional vignettes:

The rSPDE package also has a separate group of functions for performing the operator-based rational approximations introduced in Bolin and Kirchner (2020). These are especially 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), Bolin (2013) and Asar et al. (2020) 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:

A toy data set

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)
n_loc <- 2000
loc_2d_mesh <- matrix(runif(n_loc * 2), n_loc, 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 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 <- 2
range <- 0.25
nu <- 1.3
kappa <- sqrt(8 * nu) / range
op <- matern.operators(
  mesh = mesh_2d, nu = nu,
  range = range, sigma = sigma, m = 2,
  parameterization = "matern"
)
tau <- op$tau

We can now use the simulate function to simulate a realization of the field \(u\):

u <- simulate(op)

Let us now consider a simple Gaussian linear model where the spatial field \(u(\mathbf{s})\) is observed at \(m\) locations, \(\{\mathbf{s}_1 , \ldots , \mathbf{s}_m \}\) under Gaussian measurement noise. 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.

To generate a data set y from this model, we first draw some observation locations at random in the domain and then use the R-INLA function inla.spde.make.A() to construct the observation matrix which can be used to evaluate the simulated field \(u\) at the observation locations. After this we simply add the measurment noise.

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

The generated data can be seen in the following image.

library(ggplot2)
library(viridis)
#> Loading required package: viridisLite
df <- data.frame(x1 = as.double(loc_2d_mesh[, 1]),
  x2 = as.double(loc_2d_mesh[, 2]), y = as.double(y))
ggplot(df, aes(x = x1, y = x2, col = y)) +
  geom_point() +
  scale_color_viridis()

The simulated random field is shown in the following figure.

proj <- inla.mesh.projector(mesh_2d, dims = c(100, 100))
field <- inla.mesh.project(proj, field = as.vector(u))
field.df <- data.frame(x1 = proj$lattice$loc[,1],
                        x2 = proj$lattice$loc[,2], 
                        y = as.vector(field))
ggplot(field.df, aes(x = x1, y = x2, fill = y)) +
  geom_raster() + xlim(0,1) + ylim(0,1) + 
  scale_fill_viridis()

Fitting the model with R-INLA implementation of the rational SPDE approach

We 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\). The default value for this is \(4\). If we increase the upper bound for \(\nu\) we also increase the computational cost, and if we decrease the upper bound we also decrease the computatoinal cost. For this example we set nu.upper.bound=2. See the R-INLA implementation of the rational SPDE approach for further details.

rspde_model <- rspde.matern(
  mesh = mesh_2d,
  nu.upper.bound = 2,
  parameterization = "spde"
)

Finally, we create the formula and fit the model to the data:

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))
  )

We can get a summary of the fit:

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, ", " keep = keep, working.directory = working.directory, 
#>    silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = 
#>    debug, .parent.frame = .parent.frame)" ) 
#> Time used:
#>     Pre = 2.38, Running = 29.3, Post = 0.183, Total = 31.8 
#> Random effects:
#>   Name     Model
#>     field CGeneric
#> 
#> Model hyperparameters:
#>                                            mean    sd 0.025quant 0.5quant
#> Precision for the Gaussian observations 101.663 3.500     94.887  101.624
#> Theta1 for field                         -4.600 0.164     -4.929   -4.598
#> Theta2 for field                          2.425 0.160      2.104    2.428
#> Theta3 for field                          0.127 0.082     -0.031    0.125
#>                                         0.975quant   mode
#> Precision for the Gaussian observations    108.665 101.59
#> Theta1 for field                            -4.283  -4.59
#> Theta2 for field                             2.731   2.44
#> Theta3 for field                             0.292   0.12
#> 
#> Marginal log-Likelihood:  854.64 
#>  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)')

To get a summary of the fit of the random field only, we can do the following:

result_fit <- rspde.result(rspde_fit, "field", rspde_model)
summary(result_fit)
#>             mean         sd 0.025quant   0.5quant 0.975quant        mode
#> tau    0.0101871 0.00166572 0.00723344  0.0100549  0.0137771  0.00988217
#> kappa 11.4453000 1.81543000 8.21709000 11.3370000 15.3284000 11.13980000
#> nu     1.0631900 0.04065690 0.98502700  1.0625000  1.1445500  1.05987000
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.004452908  0.01018715  0.009882171
#> 2     kappa 12.899612397 11.44530109 11.139796634
#> 3        nu  1.300000000  1.06319253  1.059866434

We can also obtain the summary in the matern parameterization by setting the parameterization argument to matern:

result_fit_matern <- rspde.result(rspde_fit, "field", rspde_model,
                                  parameterization = "matern")
summary(result_fit_matern)
#>             mean        sd 0.025quant 0.5quant 0.975quant     mode
#> std.dev 2.797860 0.3854900   2.130620 2.769600   3.639850 2.733570
#> range   0.261992 0.0504719   0.178158 0.257496   0.374983 0.239651
#> nu      1.063190 0.0406569   0.985027 1.062500   1.144550 1.059870
result_df_matern <- data.frame(
  parameter = c("sigma", "range", "nu"),
  true = c(sigma, range, nu), mean = c(
    result_fit_matern$summary.std.dev$mean,
    result_fit_matern$summary.range$mean,
    result_fit_matern$summary.nu$mean
  ),
  mode = c(
    result_fit_matern$summary.std.dev$mode,
    result_fit_matern$summary.range$mode,
    result_fit_matern$summary.nu$mode
  )
)
print(result_df_matern)
#>   parameter true      mean      mode
#> 1     sigma 2.00 2.7978631 2.7335662
#> 2     range 0.25 0.2619917 0.2396512
#> 3        nu 1.30 1.0631925 1.0598664

Kriging with R-INLA implementation of the rational SPDE approach

Let 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 of locations where we want to compute 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:

projgrid <- rspde.mesh.projector(mesh_2d,
  xlim = c(0, 1),
  ylim = c(0, 1)
)

This lattice contains 100 × 100 locations (the default) which are shown in the following figure:

coord.prd <- projgrid$lattice$loc
plot(coord.prd, type = "p", cex = 0.1)

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

A.prd <- projgrid$proj$A

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:

field.pred.df <- data.frame(x1 = projgrid$lattice$loc[,1],
                        x2 = projgrid$lattice$loc[,2], 
                        y = as.vector(m.prd))
ggplot(field.pred.df, aes(x = x1, y = x2, fill = y)) +
  geom_raster() + xlim(0,1) + ylim(0,1) + 
  scale_fill_viridis()
#> Warning: Removed 396 rows containing missing values (`geom_raster()`).

Then, the marginal standard deviations:

field.pred.sd.df <- data.frame(x1 = proj$lattice$loc[,1],
                        x2 = proj$lattice$loc[,2], 
                        sd = as.vector(sd.prd))
ggplot(field.pred.sd.df, aes(x = x1, y = x2, fill = sd)) +
  geom_raster() + xlim(0,1) + ylim(0,1) + 
  geom_raster() +
  scale_fill_viridis()
#> Warning: Removed 6156 rows containing missing values (`geom_raster()`).
#> Removed 6156 rows containing missing values (`geom_raster()`).

Fitting the model with inlabru implementation of the rational SPDE approach

We will now fit the same model of the toy data set using our inlabru implementation of the rational SPDE approach. Further details on this implementation can be found in inlabru implementation of the rational SPDE approach.

We begin by loading the inlabru package:

library(inlabru)

The creation of the model object is the same as in R-INLA’s case:

rspde_model <- rspde.matern(
  mesh = mesh_2d,
  nu.upper.bound = 2,
  parameterization = "spde"
)

The advantage with inlabru is that we do not need to form the stack manually, but can simply collect the required data in a data.frame():

toy_df <- data.frame(coord1 = loc_2d_mesh[,1],
                     coord2 = loc_2d_mesh[,2],
                     y = as.vector(y))
coordinates(toy_df) <- c("coord1", "coord2")

Finally, we create the component and fit:

cmp <-
  y ~ -1 + field(coordinates, model = rspde_model)

rspde_bru_fit <-
  bru(cmp,
      data=toy_df,
    options=list(
    family = "gaussian")
  )

At this stage, we can get a summary of the fit just as in the R-INLA case:

summary(rspde_bru_fit)
#> inlabru version: 2.8.0
#> INLA version: 23.06.25
#> Components:
#> field: main = cgeneric(coordinates), group = exchangeable(1L), replicate = iid(1L)
#> Likelihoods:
#>   Family: 'gaussian'
#>     Data class: 'SpatialPointsDataFrame'
#>     Predictor: y ~ .
#> Time used:
#>     Pre = 1.96, Running = 23.1, Post = 0.599, Total = 25.6 
#> Random effects:
#>   Name     Model
#>     field CGeneric
#> 
#> Model hyperparameters:
#>                                            mean    sd 0.025quant 0.5quant
#> Precision for the Gaussian observations 101.656 3.506     94.844  101.625
#> Theta1 for field                         -4.273 0.595     -5.454   -4.270
#> Theta2 for field                          2.358 0.182      1.987    2.363
#> Theta3 for field                         -0.056 0.328     -0.696   -0.058
#>                                         0.975quant    mode
#> Precision for the Gaussian observations    108.644 101.626
#> Theta1 for field                            -3.111  -4.256
#> Theta2 for field                             2.702   2.383
#> Theta3 for field                             0.596  -0.066
#> 
#> Deviance Information Criterion (DIC) ...............: -3266.72
#> Deviance Information Criterion (DIC, saturated) ....: 2301.28
#> Effective number of parameters .....................: 297.97
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: -3258.62
#> Effective number of parameters .................: 265.28
#> 
#> Marginal log-Likelihood:  856.36 
#>  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 also obtain a summary of the field only:

result_fit <- rspde.result(rspde_bru_fit, "field", rspde_model)
summary(result_fit)
#>            mean       sd 0.025quant   0.5quant 0.975quant        mode
#> tau    0.016576 0.010492 0.00428241  0.0139956  0.0441599  0.00982083
#> kappa 10.746800 1.931920 7.31803000 10.6350000 14.8736000 10.45940000
#> nu     0.972803 0.159001 0.66779100  0.9709750  1.2868100  0.96530500
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.004452908  0.01657597  0.009820831
#> 2     kappa 12.899612397 10.74684193 10.459355647
#> 3        nu  1.300000000  0.97280349  0.965305493

Let us obtain a summary in the matern parameterization by setting the parameterization argument to matern:

result_fit_matern <- rspde.result(rspde_bru_fit, "field", rspde_model,
                                  parameterization = "matern")
summary(result_fit_matern)
#>             mean        sd 0.025quant 0.5quant 0.975quant     mode
#> std.dev 2.708580 0.3599210   2.106710 2.670470   3.531620 2.623230
#> range   0.271230 0.0426567   0.202809 0.265431   0.371138 0.252573
#> nu      0.972803 0.1590010   0.667791 0.970975   1.286810 0.965305
result_df_matern <- data.frame(
  parameter = c("sigma", "range", "nu"),
  true = c(sigma, range, nu), mean = c(
    result_fit_matern$summary.std.dev$mean,
    result_fit_matern$summary.range$mean,
    result_fit_matern$summary.nu$mean
  ),
  mode = c(
    result_fit_matern$summary.std.dev$mode,
    result_fit_matern$summary.range$mode,
    result_fit_matern$summary.nu$mode
  )
)
print(result_df_matern)
#>   parameter true      mean      mode
#> 1     sigma 2.00 2.7085799 2.6232345
#> 2     range 0.25 0.2712303 0.2525731
#> 3        nu 1.30 0.9728035 0.9653055

Kriging with inlabru implementation of the rational SPDE approach

Let 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 of the locations where we want to evaluate the predictions. We begin by creating a regular grid in and then extract the coorinates:

pred_coords <- data.frame(x1 = projgrid$lattice$loc[,1],
                          x2 = projgrid$lattice$loc[,2])
coordinates(pred_coords) <- c("x1", "x2")

Let us now compute the predictions. An advantage with inlabru is that we can do this after fitting the model to the data:

field_pred <- predict(rspde_bru_fit, pred_coords, ~field)

The following figure shows the mean of these predictions:

field_pred_df <- field_pred@data
field_pred_df <- cbind(field_pred_df, field_pred@coords)
ggplot(field_pred_df, aes(x = x1, y = x2, fill = mean)) +
  geom_raster() + xlim(0,1) + ylim(0,1) + 
  scale_fill_viridis()

The following figure shows the marginal standard deviations of the predictions:

ggplot(field_pred_df, aes(x = x1, y = x2, fill = sd)) +
  geom_raster() + xlim(0,1) + ylim(0,1) + 
  scale_fill_viridis()

An alternative and very simple approach is to use the pixels() function:

pxl <- pixels(mesh_2d)

field_pred <- predict(rspde_bru_fit, pxl, ~field)

ggplot() + gg(field_pred) +
  scale_fill_viridis() + xlim(0,1) + ylim(0,1)

Fitting the model with rSPDE

We will now fit the model of the toy data set without using R-INLA or inlabru. 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 use the function rSPDE.construct.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.

Let us create an object for estimation, a data.frame with the data and then fit the model using the rspde_lme() function.

op_est <- matern.operators(
  mesh = mesh_2d, m = 2
)

toy_df_rspde <- data.frame(coord1 = loc_2d_mesh[,1],
                     coord2 = loc_2d_mesh[,2],
                     y = as.vector(y))
fit_rspde <- rspde_lme(y ~ -1, data = toy_df_rspde, loc = c("coord1", "coord2"),
                      model = op_est, parallel = TRUE)

We can obtain the summary:

summary(fit_rspde)
#> 
#> Latent model - Whittle-Matern
#> 
#> Call:
#> rspde_lme(formula = y ~ -1, loc = c("coord1", "coord2"), data = toy_df_rspde, 
#>     model = op_est, parallel = TRUE)
#> 
#> No fixed effects.
#> 
#> Random effects:
#>        Estimate Std.error z-value
#> alpha 2.267e+00 3.858e-02  58.767
#> tau   4.663e-03 7.155e-04   6.518
#> kappa 1.347e+01 1.487e+00   9.057
#> 
#> Random effects (Matern parameterization):
#>       Estimate Std.error z-value
#> nu     1.26723   0.03858  32.847
#> sigma  1.99073   0.18256  10.905
#> range  0.23634   0.02559   9.234
#> 
#> Measurement error:
#>          Estimate Std.error z-value
#> std. dev  0.09921   0.00170   58.35
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
#> 
#> Log-Likelihood:  872.7136 
#> Number of function calls by 'optim' = 77
#> Optimization method used in 'optim' = L-BFGS-B
#> 
#> Time used to:     fit the model =  35.92145 secs 
#>   set up the parallelization = 1.77779 secs

Let us compare with the true values:

print(data.frame(
  sigma = c(sigma, fit_rspde$matern_coeff$random_effects[2]), 
  range = c(range, fit_rspde$matern_coeff$random_effects[3]),
  nu = c(nu, fit_rspde$matern_coeff$random_effects[1]),
  row.names = c("Truth", "Estimates")
))
#>              sigma     range       nu
#> Truth     2.000000 0.2500000 1.300000
#> Estimates 1.990734 0.2363381 1.267225

# Time to fit
print(fit_rspde$fitting_tim)
#> Time difference of 35.92146 secs

Kriging with rSPDE

We will now do kriging 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 data.frame with the prediction locations:

predgrid <- inla.mesh.projector(mesh_2d,
  xlim = c(0, 1),
  ylim = c(0, 1)
)
pred_coords <- data.frame(coord1 = predgrid$lattice$loc[,1],
                          coord2 = predgrid$lattice$loc[,2])

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.

pred.rspde <- predict(fit_rspde,
  data = pred_coords, loc = c("coord1", "coord2"),
  compute_variances = TRUE
)

Finally, we plot the results. First the mean:

field.pred2.df <- data.frame(x1 = predgrid$lattice$loc[,1],
                             x2 = predgrid$lattice$loc[,2],
                             y = as.vector(pred.rspde$mean))
ggplot(field.pred2.df, aes(x = x1, y = x2, fill = y)) +
  geom_raster() + xlim(0,1) + ylim(0,1) + 
  scale_fill_viridis()
#> Warning: Removed 396 rows containing missing values (`geom_raster()`).

Then, the standard deviations:

field.pred2.sd.df <-field.pred2.df <- data.frame(x1 = predgrid$lattice$loc[,1],
                             x2 = predgrid$lattice$loc[,2],
                             sd = as.vector(sqrt(pred.rspde$variance)))
ggplot(field.pred2.sd.df, aes(x = x1, y = x2, fill = sd)) +
  geom_raster() +
  scale_fill_viridis()

References

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.
Bolin, David, Alexandre B. Simas, and Zhen Xiong. 2023. “Covariance-Based Rational Approximations of Fractional SPDEs for Computationally Efficient Bayesian Inference.” Journal of Computational and Graphical Statistics.
Wallin, Jonas, and David Bolin. 2015. “Geostatistical Modelling Using Non-Gaussian Matérn Fields.” Scandinavian Journal of Statistics 42 (3): 872–90.