Skip to contents

The rSPDE package

CRAN_Status_Badge CRAN_Downloads R-CMD-check

rSPDE is an R package used for computing rational approximations of fractional SPDEs. These rational approximations can be used for computatially efficient statistical inference.

Basic statistical operations such as likelihood evaluations and kriging predictions using the fractional approximations are also implemented. The package also contains an interface to R-INLA.

Introduction

Several popular Gaussian random field models can be represented as solutions to stochastic partial differential equations (SPDEs) of the form

Lβ(τu)=𝒲. L^{\beta} (\tau u) = \mathcal{W}.

Here 𝒲\mathcal{W} is a Gaussian white noise, LL is a second-order differential operator, the fractional power β\beta determines the smoothness of uu, and τ\tau scales the variance of uu. The simplest example is a model on ℝd\mathbb{R}^d with L=κ2−ΔL = \kappa^2 - \Delta, which results in a Gaussian random field uu with a Matérn covariance function

C(h)=σ22ν−1Γ(ν)(κh)νKν(κh). C(h) = \dfrac{ \sigma^2 }{ 2 ^ {\nu-1} \Gamma (\nu)} (\kappa h)^{\nu} K_{\nu} (\kappa h).

If 2β2 \beta is an integer and if the domain 𝒟\mathcal{D} where the model is defined is bounded, then uu can be approximated by a Gaussian Markov random field (GMRF) u\mathbf { \mathrm{u} } via a finite element method (FEM) for the SPDE. Specifically, the approximation can be written as

uh(s)=∑i=1nuiφi(s). u_h(s) = \sum_{i=1}^n u_i \varphi_i (s).

Here {φi}\{\varphi_i\} are piecewise linear basis functions defined by some triangulation of 𝒟\mathcal{D} and the vector of weights u=(u1,…,un)⊤\mathbf{\mathrm{u}} = (u_1,\ldots,u_n)^\top is normally distributed, N(u,Q̃−1)N(\mathbf{\mathrm{u}}, \tilde{\mathbf{\mathrm{Q}}}^{-1}), where Q̃\tilde{ \mathbf{ \mathrm{Q} } } is sparse. See An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach for further details.

The rSPDE package provides corresponding computationally efficient approximations for the case when β\beta is a general fractional power. The main idea is to combine the FEM approximation with a rational approximation of the fractional operator. As a result, one can easily do inference and prediction using fractional SPDE models such as

(κ2−Δ)βu=𝒲. (\kappa^2-\Delta)^\beta u = \mathcal{W}.

In particular, it allows for bayesian inference of all model parameters, including the fractional parameter β\beta.

For illustration purposes, the package contains a simple FEM implementation for models on R. See the An introduction to the rSPDE package vignette for an introduction to the package. The Rational approximation with the rSPDE package and Operator-based rational approximation vignettes provide introductions to how to create and fit rSPDE models. For an introduction to the R-INLA implementation of the rSPDE models see the R-INLA implementation of the rational SPDE approach. The rSPDE documentation contains descriptions and examples of the functions in the rSPDE package.

Installation instructions

The latest CRAN release of the package can be installed directly from CRAN with install.packages("rSPDE"). The latest stable version (which is sometimes slightly more recent than the CRAN version), can be installed by using the command

remotes::install_github("davidbolin/rspde", ref = "stable")

in R. The development version can be installed using the command

remotes::install_github("davidbolin/rspde", ref = "devel")

If you want to install the package using the {r}emotes::install_github-method on Windows, you first need to install Rtools and add the paths to Rtools and gcc to the Windows PATH environment variable. This can be done for the current R session only using the commands

rtools = "C:\Rtools\bin"
gcc = "C:\Rtools\gcc-4.6.3\bin"
Sys.setenv(PATH = paste(c(gcc, rtools, Sys.getenv("PATH")), collapse = ";"))

where the variables {r}tools and gcc need to be changed if Rtools is not installed directly on C:, and gcc’s version might need to be changed depending on the version of Rtools.

Example with rSPDE - INLA

We will illustrate the rSPDE package with a kriging example using our R-INLA interface to rSPDE.

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 not analyse the full spatio-temporal data set, but instead look at the total precipitation in January

For further details on the dataset and on the commands we refer the reader to the rSPDE-INLA Vignette.

library(rSPDE)
library(ggplot2)
library(INLA)
library(splancs)
library(viridis)

#Load the data
data(PRprec)
data(PRborder)

#Get the precipitation in January
Y <- rowMeans(PRprec[, 3 + 1:31])

#Treat the data and plot
ind <- !is.na(Y)
Y <- Y[ind]
coords <- as.matrix(PRprec[ind, 1:2])
alt <- PRprec$Altitude[ind]

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

#Get distance from the sea
seaDist <- apply(spDists(coords, PRborder[1034:1078, ], longlat = TRUE), 1, 
                 min)
                 
#Create the mesh
library(fmesher)
prdomain <- fm_nonconvex_hull(coords, -0.03, -0.05, resolution = c(100, 100))
prmesh <- fm_mesh_2d(boundary = prdomain, max.edge = c(0.45, 1), cutoff = 0.2)
plot(prmesh, asp = 1, main = "")
lines(PRborder, col = 3)
points(coords[, 1], coords[, 2], pch = 19, cex = 0.5, col = "red")

#Create the observation matrix
Abar <- rspde.make.A(mesh = prmesh, loc = coords)

#Create the rspde model object
rspde_model <- rspde.matern(mesh = prmesh)

#Create the index and inla.stack object
mesh.index <- rspde.make.index(name = "field", mesh = prmesh)
stk.dat <- inla.stack(
  data = list(y = Y), A = list(Abar, 1), tag = "est",
  effects = list(
    c(
      mesh.index
    ),
    list(
      seaDist = inla.group(seaDist),
      Intercept = 1
    )
  )
)
                      
#Create the formula object and fit the model
f.s <- y ~ -1 + Intercept + f(seaDist, model = "rw1") +
  f(field, model = rspde_model)
  
rspde_fit <- inla(f.s, 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))
            
summary(rspde_fit)
## Time used:
##     Pre = 0.42, Running = 2.89, Post = 0.0521, Total = 3.36 
## Fixed effects:
##            mean   sd 0.025quant 0.5quant 0.975quant  mode kld
## Intercept 1.941 0.04      1.862    1.941       2.02 1.941   0
## 
## Random effects:
##   Name     Model
##     seaDist RW1 model
##    field CGeneric
## 
## Model hyperparameters:
##                                                    mean      sd 0.025quant
## Precision-parameter for the Gamma observations   14.682    1.05      12.72
## Precision for seaDist                          7795.411 4497.75    2422.19
## Theta1 for field                                 -0.905    4.71     -10.07
## Theta2 for field                                  1.473    1.06      -0.64
## Theta3 for field                                 -0.517    4.08      -8.65
##                                                0.5quant 0.975quant     mode
## Precision-parameter for the Gamma observations   14.647      16.84   14.585
## Precision for seaDist                          6735.175   19461.60 5066.420
## Theta1 for field                                 -0.943       8.49   -1.103
## Theta2 for field                                  1.481       3.54    1.514
## Theta3 for field                                 -0.484       7.43   -0.346
## 
## Marginal log-Likelihood:  -1252.89 
##  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)')
result_fit <- rspde.result(rspde_fit, "field", rspde_model)
summary(result_fit)
##              mean          sd  0.025quant 0.5quant 0.975quant        mode
## tau   2992.260000 5.44607e+04 9.51685e-05 0.386598 4567.78000 9.46756e-08
## kappa    7.517610 9.80636e+00 5.37556e-01 4.404580   33.87380 1.37896e+00
## nu       0.913089 7.96933e-01 3.86897e-04 0.766163    1.99873 1.35147e-06
#Plot the posterior densities
posterior_df_fit <- gg_df(result_fit)

ggplot(posterior_df_fit) + geom_line(aes(x = x, y = y)) + 
facet_wrap(~parameter, scales = "free") + labs(y = "Density")

#Create a grid to predict
nxy <- c(150, 100)
projgrid <- rspde.mesh.projector(prmesh,
  xlim = range(PRborder[, 1]),
  ylim = range(PRborder[, 2]), dims = nxy
)
xy.in <- inout(projgrid$lattice$loc, cbind(PRborder[, 1], PRborder[, 2]))
coord.prd <- projgrid$lattice$loc[xy.in, ]

#Compute A matrix and seaDist at predict locations and build the stack
A.prd <- projgrid$proj$A[xy.in, ]
seaDist.prd <- apply(spDists(coord.prd, 
    PRborder[1034:1078, ], longlat = TRUE), 1, min)
ef.prd = list(c(mesh.index), 
    list(long = inla.group(coord.prd[, 
    1]), lat = inla.group(coord.prd[, 2]), 
    seaDist = inla.group(seaDist.prd),
    Intercept = 1))
stk.prd <- inla.stack(data = list(y = NA), 
    A = list(A.prd, 1), tag = "prd", 
    effects = ef.prd)
stk.all <- inla.stack(stk.dat, stk.prd)

rspde_fitprd <- inla(f.s,
  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")
)

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

#Plot the predictions
pred_df <- data.frame(x1 = coord.prd[,1],
                      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()

Example with rSPDE - inlabru

We will now illustrate the rSPDE the same kriging example above with our inlabru interface to rSPDE. We will make this description self-contained, so we will not use any information or codes from the example above.

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 not analyse the full spatio-temporal data set, but instead look at the total precipitation in January

For further details on the dataset and on the commands we refer the reader to the rSPDE-inlabru Vignette.

library(rSPDE)
library(ggplot2)
library(INLA)
library(inlabru)
library(splancs)

#Load the data
data(PRprec)
data(PRborder)

#Get the precipitation in January
Y <- rowMeans(PRprec[, 3 + 1:31])

#Treat the data and plot
ind <- !is.na(Y)
Y <- Y[ind]
coords <- as.matrix(PRprec[ind, 1:2])
alt <- PRprec$Altitude[ind]

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

#Get distance from the sea
seaDist <- apply(spDists(coords, PRborder[1034:1078, ], longlat = TRUE), 1, 
                 min)
                 
#Create the mesh
library(fmesher)
prdomain <- fm_nonconvex_hull(coords, -0.03, -0.05, resolution = c(100, 100))
prmesh <- fm_mesh_2d(boundary = prdomain, max.edge = c(0.45, 1), cutoff = 0.2)
plot(prmesh, asp = 1, main = "")
lines(PRborder, col = 3)
points(coords[, 1], coords[, 2], pch = 19, cex = 0.5, col = "red")

## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
#Create the rspde model object
rspde_model <- rspde.matern(mesh = prmesh)

#Create the data.frame
prdata <- data.frame(long = coords[,1], lat = coords[,2], 
                        seaDist = inla.group(seaDist), y = Y)
                      
prdata <- st_as_sf(prdata, coords = c("long", "lat"), crs = 4326)                      
#Create the component

cmp <- y ~ Intercept(1) + distSea(seaDist, model="rw1") +
field(geometry, model = rspde_model)

# Fit the model
  
rspde_fit <- bru(cmp, family = "Gamma", 
            data = prdata,
            options = list(
            verbose = FALSE, 
            control.inla=list(int.strategy='eb'),
            control.predictor = list(compute = TRUE))
)
            
summary(rspde_fit)
## inlabru version: 2.11.1
## INLA version: 24.11.07-4
## Components:
## Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
## distSea: main = rw1(seaDist), group = exchangeable(1L), replicate = iid(1L)
## field: main = cgeneric(geometry), group = exchangeable(1L), replicate = iid(1L)
## Likelihoods:
##   Family: 'Gamma'
##     Data class: 'sf', 'data.frame'
##     Predictor: y ~ .
## Time used:
##     Pre = 0.236, Running = 3.58, Post = 0.0629, Total = 3.88 
## Fixed effects:
##            mean   sd 0.025quant 0.5quant 0.975quant  mode kld
## Intercept 1.941 0.04      1.862    1.941       2.02 1.941   0
## 
## Random effects:
##   Name     Model
##     distSea RW1 model
##    field CGeneric
## 
## Model hyperparameters:
##                                                   mean       sd 0.025quant
## Precision-parameter for the Gamma observations   14.69    1.049      12.72
## Precision for distSea                          7749.25 4463.693    2386.05
## Theta1 for field                                 -4.88    3.896     -14.01
## Theta2 for field                                  2.20    0.793       1.01
## Theta3 for field                                  2.88    3.347      -1.86
##                                                0.5quant 0.975quant     mode
## Precision-parameter for the Gamma observations    14.65   1.69e+01   14.588
## Precision for distSea                           6704.11   1.93e+04 5042.968
## Theta1 for field                                  -4.31   6.36e-01   -1.315
## Theta2 for field                                   2.09   4.03e+00    1.566
## Theta3 for field                                   2.40   1.07e+01   -0.151
## 
## Deviance Information Criterion (DIC) ...............: 2466.03
## Deviance Information Criterion (DIC, saturated) ....: 723.98
## Effective number of parameters .....................: 111.27
## 
## Watanabe-Akaike information criterion (WAIC) ...: 2467.24
## Effective number of parameters .................: 97.00
## 
## Marginal log-Likelihood:  -1253.44 
##  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)')
#Get the summary on the user's scale
result_fit <- rspde.result(rspde_fit, "field", rspde_model)
## Warning in rspde.result(rspde_fit, "field", rspde_model): the mean or mode of
## nu is very close to nu.upper.bound, please consider increasing nu.upper.bound,
## and refitting the model.
summary(result_fit)
##            mean        sd  0.025quant  0.5quant 0.975quant        mode
## tau    0.239787  0.684428 4.88542e-10 0.0152431    1.91226 4.88542e-10
## kappa 13.007200 15.976000 2.77104e+00 7.9010500   55.13660 4.07965e+00
## nu     1.523390  0.563400 2.64536e-01 1.8149100    1.99989 1.99999e+00
#Plot the posterior densities
posterior_df_fit <- gg_df(result_fit)

ggplot(posterior_df_fit) + geom_line(aes(x = x, y = y)) + 
facet_wrap(~parameter, scales = "free") + labs(y = "Density")

#Create a grid to predict
nxy <- c(150, 100)
projgrid <- rspde.mesh.projector(prmesh, xlim = range(PRborder[, 1]), 
ylim = range(PRborder[,2]), dims = nxy)
xy.in <- inout(projgrid$lattice$loc, cbind(PRborder[, 1], PRborder[, 2]))
coord.prd <- projgrid$lattice$loc[xy.in, ]

#Compute seaDist at predict locations 
seaDist.prd <- apply(spDists(coord.prd, 
    PRborder[1034:1078, ], longlat = TRUE), 1, min)

# Build the prediction data.frame()
coord.prd.df <- data.frame(long = coord.prd[,1],
                            lat = coord.prd[,2])

coord.prd.df$seaDist <- seaDist.prd

coord.prd.df <- st_as_sf(coord.prd.df, coords = c("long", "lat"), 
                        crs = 4326)

# Obtain prediction at the locations
pred_obs <- predict(rspde_fit, coord.prd.df, 
        ~exp(Intercept + field + distSea))

Finally, we plot the results. First the predicted mean:

ggplot() + gg(pred_obs, geom = "tile", 
          aes(fill = mean)) +
  scale_fill_viridis()

Then, the std. deviations:

ggplot() + gg(pred_obs, geom = "tile", aes(fill = sd)) +
  geom_raster() + scale_fill_viridis()