## Setting things up

``````library(INLA)
library(inlabru)
library(fmesher)
library(mgcv)
library(ggplot2)
library(sp)``````

Make a shortcut to a nicer colour scale:

``````colsc <- function(...) {
colours = rev(RColorBrewer::brewer.pal(11, "RdYlBu")),
limits = range(..., na.rm = TRUE)
)
}``````

## Modelling on 2D domains

We will now construct a 2D model, generate a sample of a random field, and attempt to recover the field from observations at a few locations. Tomorrow, we will look into more general mesh constructions that adapt to irregular domains.

First, we build a high resolution mesh for the true field, using low level INLA functions

``````bnd <- spoly(data.frame(easting = c(0, 10, 10, 0), northing = c(0, 0, 10, 10)))
mesh_fine <- fm_mesh_2d_inla(boundary = bnd, max.edge = 0.2)
ggplot() +
geom_fm(data = mesh_fine)
#> Warning in fm_as_sfc.fm_segm(data): fm_as_sfc currently only supports
#> (multi)linestring output``````

``````
# Note: the priors here will not be used in estimation
matern_fine <-
inla.spde2.pcmatern(mesh_fine,
prior.sigma = c(1, 0.01),
prior.range = c(1, 0.01)
)
true_range <- 4
true_sigma <- 1
true_Q <- inla.spde.precision(matern_fine, theta = log(c(true_range, true_sigma)))``````

What is the pointwise standard deviation of the field? Along straight boundaries, the variance is twice the target variance. At corners the variance is 4 times as large.

``````true_sd <- diag(inla.qinv(true_Q))^0.5
ggplot() +
gg(mesh_fine, color = true_sd) +
coord_equal()``````

Generate a sample from the model:

``````true_field <- inla.qsample(1, true_Q)[, 1]

truth <- expand.grid(
easting = seq(0, 10, length = 100),
northing = seq(0, 10, length = 100)
)
truth <- sf::st_as_sf(truth, coords = c("easting", "northing"))
truth\$field <- fm_evaluate(
mesh_fine,
loc = truth,
field = true_field
)

pl_truth <- ggplot() +
gg(truth, aes(fill = field), geom = "tile") +
ggtitle("True field")
pl_truth``````

``````
## Or with another colour scale:
csc <- colsc(truth\$field)
multiplot(pl_truth, pl_truth + csc, cols = 2)``````

Extract observations from some random locations:

``````n <- 200
mydata <- sf::st_as_sf(
data.frame(easting = runif(n, 0, 10), northing = runif(n, 0, 10)),
coords = c("easting", "northing")
)
mydata\$observed <-
fm_evaluate(
mesh_fine,
loc = mydata,
field = true_field
) +
rnorm(n, sd = 0.4)
ggplot() +
gg(mydata, aes(col = observed))``````

## Estimating the field

Construct a mesh covering the data:

``````mesh <- fm_mesh_2d_inla(boundary = bnd, max.edge = 0.5)
ggplot() +
geom_fm(data = mesh)``````

Construct an SPDE model object for a Matern model:

``````matern <-
inla.spde2.pcmatern(mesh,
prior.sigma = c(10, 0.01),
prior.range = c(1, 0.01)
)``````

Specify the model components:

``cmp <- observed ~ field(geometry, model = matern) + Intercept(1)``

Fit the model and inspect the results:

``````fit <- bru(cmp, mydata, family = "gaussian")
summary(fit)``````

Predict the field on a lattice, and generate a single realisation from the posterior distribution:

``````pix <- fm_pixels(mesh, dims = c(200, 200))
pred <- predict(
fit, pix,
~ field + Intercept
)
samp <- generate(fit, pix,
~ field + Intercept,
n.samples = 1
)
pred\$sample <- samp[, 1]``````

Compare the truth to the estimated field (posterior mean and a sample from the posterior distribution):

``````pl_posterior_mean <- ggplot() +
gg(pred, geom = "tile") +
gg(bnd, alpha = 0) +
ggtitle("Posterior mean")
pl_posterior_sample <- ggplot() +
gg(pred, aes(fill = sample), geom = "tile") +
gg(bnd, alpha = 0) +
ggtitle("Posterior sample")

# Common colour scale for the truth and estimate:
csc <- colsc(truth\$field, pred\$mean, pred\$sample)
multiplot(pl_truth + csc,
pl_posterior_mean + csc,
pl_posterior_sample + csc,
cols = 3
)``````

Plot the SPDE parameter and fixed effect parameter posteriors.

``````int.plot <- plot(fit, "Intercept")
spde.range <- spde.posterior(fit, "field", what = "range")
spde.logvar <- spde.posterior(fit, "field", what = "log.variance")
range.plot <- plot(spde.range)
var.plot <- plot(spde.logvar)

multiplot(range.plot, var.plot, int.plot)``````

Look at the correlation function if you want to:

``````corplot <- plot(spde.posterior(fit, "field", what = "matern.correlation"))
covplot <- plot(spde.posterior(fit, "field", what = "matern.covariance"))
multiplot(covplot, corplot)``````

You can plot the median, lower 95% and upper 95% density surfaces as follows (assuming that the predicted intensity is in object `pred`).

``````csc <- colsc(
pred[["median"]],
pred[["q0.025"]],
pred[["q0.975"]]
) ## Common colour scale from SpatialPixelsDataFrame

gmedian <- ggplot() +
gg(pred["median"], geom = "tile") +
csc
glower95 <- ggplot() +
gg(pred["q0.025"], geom = "tile") +
csc +
theme(legend.position = "none")
gupper95 <- ggplot() +
gg(pred["q0.975"], geom = "tile") +
csc +
theme(legend.position = "none")

multiplot(gmedian, glower95, gupper95,
layout = matrix(c(1, 1, 2, 3), byrow = TRUE, ncol = 2)
)``````