Skip to contents

Intrinsic models in the rSPDE package

David Bolin

2024-12-20

Source: vignettes/intrinsic.Rmd
intrinsic.Rmd

Introduction

In this vignette we provide a brief introduction to the intrinsic models implemented in the rSPDE package.

A fractional intrinsic model

A basic intrinsic model which is implemented in rSPDE is defined as
(−Δ)β/2(τu)=𝒲, (-\Delta)^{\beta/2}(\tau u) = \mathcal{W}, where β>d/2\beta > d/2 and dd is the dimension of the spatial domain.

To illustrate these models, we begin by defining a mesh over [0,2]×[0,2][0,2]\times [0, 2]:

library(fmesher)
bnd <- fm_segm(rbind(c(0, 0), c(2, 0), c(2, 2), c(0, 2)), is.bnd = TRUE)
mesh_2d <- fm_mesh_2d(
    boundary = bnd, 
    cutoff = 0.04,
    max.edge = c(0.1)
)
plot(mesh_2d, main = "")

We now use the intrinsic.operators() function to construct the rSPDE representation of the general model.

library(rSPDE)
#> Error in get(paste0(generic, ".", class), envir = get_method_env()) : 
#>   object 'type_sum.accel' not found
tau <- 0.2
beta <- 1.8
rspde.order <- 2
op <- intrinsic.operators(tau = tau, beta = beta, mesh = mesh_2d, m = rspde.order)

To see that the rSPDE model is approximating the true model, we can compare the variogram of the approximation with the true variogram (implemented in variogram.intrinsic.spde()) as follows.

ind.fix <- 1 + seq(from=0,to=op$m_beta * mesh_2d$n, by = mesh_2d$n)
Sigma <- op$A[,-ind.fix] %*% solve(op$Q[-ind.fix,-ind.fix], t(op$A[,-ind.fix]))
One <- rep(1, times = ncol(Sigma))
D <- diag(Sigma)
Gamma <- 0.5 * (One %*% t(D) + D %*% t(One) - 2 * Sigma)
point <- matrix(c(1,1),1,2)
Aobs <- spde.make.A(mesh = mesh_2d, loc = point)
vario <- variogram.intrinsic.spde(point, mesh_2d$loc[,1:2], tau = tau,
                                  beta = beta, L = 2, d = 2)

d = sqrt((mesh_2d$loc[,1]-point[1])^2 +  (mesh_2d$loc[,2]-point[2])^2)
plot(d, Aobs%*%Gamma, xlim = c(0,0.7))
lines(sort(d),sort(vario),col=2, lwd = 2)

We can note that the approximation is not great here, and the reason is that the order of the rational approximation is too low. Let us increase it and recompute:

rspde.order <- 8
op <- intrinsic.operators(tau = tau, beta = beta, mesh = mesh_2d, m = rspde.order)
ind.fix <- 1 + seq(from=0,to=op$m_beta * mesh_2d$n, by = mesh_2d$n)
Sigma <- op$A[,-ind.fix] %*% solve(op$Q[-ind.fix,-ind.fix], t(op$A[,-ind.fix]))
One <- rep(1, times = ncol(Sigma))
D <- diag(Sigma)
Gamma <- 0.5 * (One %*% t(D) + D %*% t(One) - 2 * Sigma)
plot(d, Aobs%*%Gamma, xlim = c(0,0.7))
lines(sort(d),sort(vario),col=2, lwd = 2)

We can now use the simulate function to simulate a realization of the field uu:

u <- simulate(op,nsim = 1)

proj <- fm_evaluator(mesh_2d, dims = c(100, 100))
field <- fm_evaluate(proj, field = as.vector(u))
field.df <- data.frame(x1 = proj$lattice$loc[,1],
                       x2 = proj$lattice$loc[,2], 
                       y = as.vector(field))

library(ggplot2)
library(viridis)
#> Loading required package: viridisLite
ggplot(field.df, aes(x = x1, y = x2, fill = y)) +
    geom_raster() +
    scale_fill_viridis()

By default, the field is simulated with a zero-integral constraint.

Fitting the model with R-INLA

Let us now consider a simple Gaussian linear model where the spatial field u(𝐬)u(\mathbf{s}) is observed at mm locations, {𝐬1,…,𝐬m}\{\mathbf{s}_1 , \ldots , \mathbf{s}_m \} under Gaussian measurement noise. For each i=1,…,m,i = 1,\ldots,m, we have yi=u(𝐬i)+εi, \begin{align} y_i &= u(\mathbf{s}_i)+\varepsilon_i\\ \end{align}, where ε1,…,εm\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 spde.make.A() functions (that wraps the functions fm_basis(), fm_block() and fm_row_kron() of the fmesher package) to construct the observation matrix which can be used to evaluate the simulated field uu at the observation locations. After this we simply add the measurment noise.

n_loc <- 1000
loc_2d_mesh <- matrix(2*runif(n_loc * 2), n_loc, 2)

A <- 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.

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

We will now fit the model 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.

library(INLA)
#> This is INLA_24.12.11 built 2024-12-11 19:58:26 UTC.
#>  - See www.r-inla.org/contact-us for how to get help.
#>  - List available models/likelihoods/etc with inla.list.models()
#>  - Use inla.doc(<NAME>) to access documentation
mesh.index <- rspde.make.index(name = "field", mesh = mesh_2d, rspde.order = rspde.order)
Abar <- rspde.make.A(mesh = mesh_2d, loc = loc_2d_mesh, rspde.order = rspde.order)
st.dat <- inla.stack(data = list(y = as.vector(y)), A = Abar, effects = mesh.index)

We now create the model object.

rspde_model <- rspde.intrinsic(mesh = mesh_2d, rspde.order = rspde.order)

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

To compare the estimated parameters to the true parameters, 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.171059 0.0227755   0.128468 0.170561   0.217571 0.170605
#> nu  0.870386 0.0371077   0.801976 0.868599   0.947309 0.862762
tau <- op$tau
nu <- op$beta - 1 #beta = nu + d/2 
result_df <- data.frame(
    parameter = c("tau", "nu", "sigma.e"),
    true = c(tau, nu, sigma.e), 
    mean = c(result_fit$summary.tau$mean,result_fit$summary.nu$mean,
             sqrt(1/rspde_fit$summary.hyperpar[1,1])),
    mode = c(result_fit$summary.tau$mode, result_fit$summary.nu$mode,
             sqrt(1/rspde_fit$summary.hyperpar[1,6]))
)
print(result_df)
#>   parameter true      mean      mode
#> 1       tau  0.2 0.1710593 0.1706048
#> 2        nu  0.8 0.8703861 0.8627616
#> 3   sigma.e  0.1 0.1073880 0.1080170

Extreme value models

When used for extreme value statistics, one might want to use a particular form of the mean value of the latent field uu, which is zero at one location kk and is given by the diagonal of Q−k,−k−1Q_{-k,-k}^{-1} for the remaining locations. This option can be specified via the mean.correction argument of rspde.intrinsic:

rspde_model2 <- rspde.intrinsic(mesh = mesh_2d, rspde.order = rspde.order,
                                mean.correction = TRUE)

We can then fit this model as before:

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

To see the posterior distributions of the parameters we can do:

result_fit <- rspde.result(rspde_fit, "field", rspde_model2)
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")

A more general model

The rSPDE package also contains a partial implementation of a more general intrinsic model, which we refer to as an intrinsic Matérn model. The model is defined as
(−Δ)β/2(κ2−Δ)α/2(τu)=𝒲, (-\Delta)^{\beta/2}(\kappa^2-\Delta)^{\alpha/2}(\tau u) = \mathcal{W}, where α+β>d/2\alpha + \beta > d/2 and dd is the dimension of the spatial domain. These models are handled by performing two rational approximations, one for each fractional operator.

To illustrate this model, we consider the same mesh as before and use the intrinsic.matern.operators() function to construct the rSPDE representation of the general model.

kappa <- 10
tau <- 0.01
alpha <- 2
beta <- 1
op <- intrinsic.matern.operators(kappa = kappa, tau = tau, alpha = alpha, 
                                 beta = beta, mesh = mesh_2d)

To see that the rSPDE model is approximating the true model, we can compare the variogram of the approximation with the true variogram (implemented in variogram.intrinsic.spde()) as follows.

Sigma <- op$A[,-1] %*% solve(op$Q[-1,-1], t(op$A[,-1]))
One <- rep(1, times = ncol(Sigma))
D <- diag(Sigma)
Gamma <- 0.5 * (One %*% t(D) + D %*% t(One) - 2 * Sigma)
point <- matrix(c(1,1),1,2)
Aobs <- spde.make.A(mesh = mesh_2d, loc = point)
vario <- variogram.intrinsic.spde(point, mesh_2d$loc[,1:2], kappa = kappa, 
                                  alpha = alpha, tau = tau,
                                  beta = beta, L = 2, d = 2)

d = sqrt((mesh_2d$loc[,1]-point[1])^2 +  (mesh_2d$loc[,2]-point[2])^2)
plot(d, Aobs%*%Gamma, xlim = c(0,0.5), ylim = c(0,0.2))
lines(sort(d),sort(vario),col=2, lwd = 2)

We can now use the simulate function to simulate a realization of the field uu:

u <- simulate(op,nsim = 1)

proj <- fm_evaluator(mesh_2d, dims = c(100, 100))
field <- fm_evaluate(proj, field = as.vector(u))
field.df <- data.frame(x1 = proj$lattice$loc[,1],
                       x2 = proj$lattice$loc[,2], 
                       y = as.vector(field))

library(ggplot2)
library(viridis)
ggplot(field.df, aes(x = x1, y = x2, fill = y)) +
    geom_raster() +
    scale_fill_viridis()

By default, the field is simulated with a zero-integral constraint.

Fitting the model with R-INLA

We will now fit the model 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 simulating some data as before.

n_loc <- 2000
loc_2d_mesh <- matrix(2*runif(n_loc * 2), n_loc, 2)

A <- 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.

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

To fit the model, we create the AA matrix, the index, and the inla.stack object. For now, these more general models can only be estimated with β=1\beta = 1 and α=1\alpha = 1 or α=2\alpha = 2. For these non-fractional models, we can use the standard INLA functions to make the required elements.

mesh.index <- inla.spde.make.index(name = "field", n.spde = mesh_2d$n)

st.dat <- inla.stack(data = list(y = as.vector(y)), A = A, effects = mesh.index)

We now create the model object.

rspde_model <- rspde.intrinsic.matern(mesh = mesh_2d, alpha = alpha)

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)
#> Time used:
#>     Pre = 0.22, Running = 3.43, Post = 0.0435, Total = 3.69 
#> Random effects:
#>   Name     Model
#>     field CGeneric
#> 
#> Model hyperparameters:
#>                                           mean    sd 0.025quant 0.5quant
#> Precision for the Gaussian observations 101.85 3.716      94.77   101.76
#> Theta1 for field                         -4.75 0.098      -4.94    -4.75
#> Theta2 for field                          2.47 0.113       2.25     2.46
#>                                         0.975quant   mode
#> Precision for the Gaussian observations     109.40 101.56
#> Theta1 for field                             -4.56  -4.74
#> Theta2 for field                              2.69   2.46
#> 
#> Marginal log-Likelihood:  1353.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)')

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.0087198 0.000847607 0.00714327  0.00867942  0.0104744  0.00863805
#> kappa 11.8560000 1.345090000 9.46727000 11.76600000 14.7439000 11.57430000
tau <- op$tau
result_df <- data.frame(
  parameter = c("tau", "kappa"),
  true = c(tau, kappa), mean = c(result_fit$summary.tau$mean,
                                     result_fit$summary.kappa$mean),
  mode = c(result_fit$summary.tau$mode, result_fit$summary.kappa$mode)
)
print(result_df)
#>   parameter  true         mean         mode
#> 1       tau  0.01  0.008719804  0.008638048
#> 2     kappa 10.00 11.855977452 11.574306686

Kriging with R-INLA implementation

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 <- inla.mesh.projector(mesh_2d,
  xlim = c(0, 2),
  ylim = c(0, 2)
)

This lattice contains 100 × 100 locations (the default). 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 AA 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()  + scale_fill_viridis()

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