An example with real data
To illustrate our implementation of rSPDE in inlabru we will consider a
dataset available in R-INLA. This data has
also been used to illustrate the SPDE approach, see for instance the
book Advanced
Spatial Modeling with Stochastic Partial Differential Equations Using R
and INLA and also the vignette Spatial
Statistics using R-INLA and Gaussian Markov random fields. See also
Lindgren, Rue, and Lindström (2011) for
theoretical details on the standard SPDE approach.
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.
An rSPDE model for precipitation
We will follow the vignette Spatial
Statistics using R-INLA and Gaussian Markov random fields. As
precipitation data are always positive, we will assume it is Gamma
distributed. R-INLA
uses the following parameterization of the Gamma distribution, \[\Gamma(\mu, \phi): \pi (y) =
\frac{1}{\Gamma(\phi)} \left(\frac{\phi}{\mu}\right)^{\phi} y^{\phi - 1}
\exp\left(-\frac{\phi y}{\mu}\right) .\] In this
parameterization, the distribution has expected value \(E(x) = \mu\) and variance \(V(x) = \mu^2/\phi\), where \(1/\phi\) is a dispersion parameter.
In this example \(\mu\) will be modelled using a stochastic model that includes both covariates and spatial structure, resulting in the latent Gaussian model for the precipitation measurements \[\begin{align} y_i\mid \mu(s_i), \theta &\sim \Gamma(\mu(s_i),c\phi)\\ \log (\mu(s)) &= \eta(s) = \sum_k f_k(c_k(s))+u(s)\\ \theta &\sim \pi(\theta) \end{align},\]
where \(y_i\) denotes the
measurement taken at location \(s_i\),
\(c_k(s)\) are covariates, \(u(s)\) is a mean-zero Gaussian Matérn
field, and \(\theta\) is a vector
containing all parameters of the model, including smoothness of the
field. That is, by using the rSPDE model we will also be
able to estimate the smoothness of the latent field.
Examining the data
We will be using inlabru. The
inlabru package is available on CRAN and also on GitHub.
We begin by loading some libraries we need to get the data and build the plots.
library(ggplot2)
library(INLA)
library(inlabru)
library(splancs)
library(viridis)Let us load the data and the border of the region
data(PRprec)
data(PRborder)The data frame contains daily measurements at 616 stations for the year 2011, as well as coordinates and altitude information for the measurement stations. We will not analyze the full spatio-temporal data set, but instead look at the total precipitation in January, which we calculate as
Y <- rowMeans(PRprec[, 3 + 1:31])In the next snippet of code, we extract the coordinates and altitudes and remove the locations with missing values.
ind <- !is.na(Y)
Y <- Y[ind]
coords <- as.matrix(PRprec[ind, 1:2])
alt <- PRprec$Altitude[ind]Let us build a plot for the precipitations:
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()
The red line in the figure shows the coast line, and we expect the distance to the coast to be a good covariate for precipitation.
This covariate is not available, so let us calculate it for each observation location:
seaDist <- apply(spDists(coords, PRborder[1034:1078, ],
longlat = TRUE
), 1, min)Now, let us plot the precipitation as a function of the possible covariates:
par(mfrow = c(2, 2))
plot(coords[, 1], Y, cex = 0.5, xlab = "Longitude")
plot(coords[, 2], Y, cex = 0.5, xlab = "Latitude")
plot(seaDist, Y, cex = 0.5, xlab = "Distance to sea")
plot(alt, Y, cex = 0.5, xlab = "Altitude")
par(mfrow = c(1, 1))Creating the rSPDE model
To use the inlabru
implementation of the rSPDE model we need to load the
functions:
library(rSPDE)To create a rSPDE model, one would the
rspde.matern() function in a similar fashion as one would
use the inla.spde2.matern() function.
Mesh
We can use R-INLA
for creating the mesh. Let us create a mesh which is based on a
non-convex hull to avoid adding many small triangles outside the domain
of interest:
prdomain <- inla.nonconvex.hull(coords, -0.03, -0.05, resolution = c(100, 100))
prmesh <- inla.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")
Setting up the data frame
In place of a inla.stack, we can set up a
data.frame() to use inlabru. We refer the reader
to vignettes in https://inlabru-org.github.io/inlabru/index.html for
further details.
prdata <- data.frame(long = coords[,1], lat = coords[,2],
seaDist = inla.group(seaDist), y = Y)
coordinates(prdata) <- c("long","lat")Setting up the rSPDE model
To set up an rSPDEmodel, all we need is the mesh. By
default it will assume that we want to estimate the smoothness parameter
\(\nu\) and to do a covariance-based
rational approximation of order 2.
Later in this vignette we will also see other options for setting up
rSPDE models such as keeping the smoothness parameter fixed
and/or increasing the order of the covariance-based rational
approximation.
Therefore, to set up a model all we have to do is use the
rspde.matern() function:
rspde_model <- rspde.matern(mesh = prmesh)Notice that this function is very reminiscent of R-INLA’s
inla.spde2.matern() function.
We will assume the following linkage between model components and observations \[\eta(s) \sim A x(s) + A \text{ Intercept} + \text{seaDist}.\] \(\eta(s)\) will then be used in the observation-likelihood, \[y_i\mid \eta(s_i),\theta \sim \Gamma(\exp(\eta (s_i)), c\phi).\]
Model fitting
We will build a model using the distance to the sea \(x_i\) as a covariate through an improper
CAR(1) model with \(\beta_{ij}=1(i\sim
j)\), which R-INLA calls a random
walk of order 1. We will fit it in inlabru’s style:
cmp <- y ~ Intercept(1) + distSea(seaDist, model="rw1") +
field(coordinates, model = rspde_model)To fit the model we simply use the bru() function:
rspde_fit <- bru(cmp, data = prdata,
family = "Gamma",
options = list(
control.inla = list(int.strategy = "eb"),
verbose = FALSE)
)inlabru results
We can look at some summaries of the posterior distributions for the parameters, for example the fixed effects (i.e. the intercept) and the hyper-parameters (i.e. dispersion in the gamma likelihood, the precision of the RW1, and the parameters of the spatial field):
summary(rspde_fit)## inlabru version: 2.7.0
## INLA version: 23.01.12
## Components:
## Intercept: main = linear(1)
## distSea: main = rw1(seaDist)
## field: main = cgeneric(coordinates)
## Likelihoods:
## Family: 'Gamma'
## Data class: 'SpatialPointsDataFrame'
## Predictor: y ~ .
## Time used:
## Pre = 3.52, Running = 11.5, Post = 0.0952, Total = 15.1
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 1.943 0.052 1.84 1.943 2.046 1.943 0
##
## Random effects:
## Name Model
## distSea RW1 model
## field CGeneric
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision parameter for the Gamma observations 13.248 0.904 11.519
## Precision for distSea 10076.588 7538.312 2450.317
## Theta1 for field -1.477 0.152 -1.787
## Theta2 for field -0.014 0.336 -0.641
## Theta3 for field -0.438 0.921 -2.146
## 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 13.232 15.08 13.226
## Precision for distSea 8012.140 30096.37 5321.865
## Theta1 for field -1.473 -1.19 -1.459
## Theta2 for field -0.027 0.68 -0.076
## Theta3 for field -0.477 1.48 -0.633
##
## Deviance Information Criterion (DIC) ...............: 2501.62
## Deviance Information Criterion (DIC, saturated) ....: 699.67
## Effective number of parameters .....................: 85.26
##
## Watanabe-Akaike information criterion (WAIC) ...: 2504.04
## Effective number of parameters .................: 77.32
##
## Marginal log-Likelihood: -1262.47
## 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)')Let \(\theta_1 = \textrm{Theta1}\), \(\theta_2=\textrm{Theta2}\) and \(\theta_3=\textrm{Theta3}\). In terms of the SPDE \[(\kappa^2 I - \Delta)^{\alpha/2}(\tau u) = \mathcal{W},\] where \(\alpha = \nu + d/2\), we have that \[\tau = \exp(\theta_1),\quad \kappa = \exp(\theta_2), \] and by default \[\nu = 4\Big(\frac{\exp(\theta_3)}{1+\exp(\theta_3)}\Big).\] The number 4 comes from the upper bound for \(\nu\), which is discussed in R-INLA implementation of the rational SPDE approach vignette.
In general, we have \[\nu = \nu_{UB}\Big(\frac{\exp(\theta_3)}{1+\exp(\theta_3)}\Big),\] where \(\nu_{UB}\) is the value of the upper bound for the smoothness parameter \(\nu\).
Another choice for prior for \(\nu\) is a truncated lognormal distribution and is also discussed in R-INLA implementation of the rational SPDE approach vignette.
inlabru results in the original scale
We can obtain outputs with respect to parameters in the original
scale by using the function rspde.result():
result_fit <- rspde.result(rspde_fit, "field", rspde_model)
summary(result_fit)## mean sd 0.025quant 0.5quant 0.975quant mode
## std.dev 0.230956 0.0348282 0.168015 0.229275 0.304398 0.226752
## range 1.043220 0.3689280 0.529886 0.971506 1.960340 0.846184
## nu 1.218190 0.5759570 0.320003 1.143500 2.434650 0.847655We can also plot the posterior densities. To this end we will use the
gg_df() function, which creates ggplot2
user-friendly data frames:
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")
Predictions
Let us now obtain predictions (i.e. do kriging) of the expected precipitation 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 inla.mesh.projector()
function:
nxy <- c(150, 100)
projgrid <- inla.mesh.projector(prmesh,
xlim = range(PRborder[, 1]),
ylim = range(PRborder[, 2]), dims = nxy
)This lattice contains 150 × 100 locations. One can easily change the
resolution of the kriging prediction by changing nxy. Let
us find the cells that are outside the region of interest so that we do
not plot the estimates there.
xy.in <- inout(projgrid$lattice$loc, cbind(PRborder[, 1], PRborder[, 2]))Let us plot the locations that we will do prediction:
coord.prd <- projgrid$lattice$loc[xy.in, ]
plot(coord.prd, type = "p", cex = 0.1)
lines(PRborder)
points(coords[, 1], coords[, 2], pch = 19, cex = 0.5, col = "red")
Let us now create a data.frame() of the coordinates:
coord.prd.df <- data.frame(x1 = coord.prd[,1],
x2 = coord.prd[,2])
coordinates(coord.prd.df) <- c("x1", "x2")Since we are using distance to the sea as a covariate, we also have
to calculate this covariate for the prediction locations. Finally, we
add the prediction location to our prediction data.frame(),
namely, coord.prd.df:
seaDist.prd <- apply(spDists(coord.prd,
PRborder[1034:1078, ],
longlat = TRUE
), 1, min)
coord.prd.df$seaDist <- seaDist.prdpred_obs <- predict(rspde_fit, coord.prd.df,
~exp(Intercept + field + distSea))Let us now build the data frame with the results:
pred_df <- pred_obs@data
pred_df <- cbind(pred_df, pred_obs@coords)Finally, we plot the results. First the predicted mean:
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()
