Example with real data
To illustrate our implementation of rSPDE in R-INLA 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)(https://rss.onlinelibrary.wiley.com/doi/full/10.1111/j.1467-9868.2011.00777.x) 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, c\phi): \pi (y) = \frac{1}{\Gamma(c\phi)} \left(\frac{c\phi}{\mu}\right)^{c\phi} y^{c\phi - 1} \exp\left(-\frac{c\phi y}{\mu}\right) .\] In this parameterization, the distribution has expected value \(E(x) = \mu\) and variance \(V(x) = \mu^2/(c\phi)\), where \(c\) is a fixed (known) scaling parameter and \(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 R-INLA. To install R-INLA go to R-INLA Project.
We begin by loading some libraries we need to get the data and build the plots.
Let us load the data and the border of the region
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
In the next snippet of code, we extract the coordinates and altitudes and remove the locations with missing values.
Let us build a plot for the precipitations:
ggplot() + geom_point(aes(x = coords[, 1], y = coords[, 2], colour = Y), size = 2,
alpha = 1) + scale_colour_gradientn(colours = tim.colors(100)) + geom_path(aes(x = PRborder[,
1], y = PRborder[, 2])) + geom_path(aes(x = PRborder[1034:1078, 1], y = PRborder[1034:1078,
2]), colour = "red")
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:
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")
Creating the rSPDE model
To use the R-INLA implementation of the rSPDE model we need to load the functions:
The rSPDE-INLA implementation is very reminiscent of R-INLA, so its usage should be straightforward for R-INLA users. For instance, to create a rSPDE model, one would rspde.matern() in place of inla.spde2.matern(). To create an index, one should use rspde.make.index() in place of inla.spde.make.index(). To create the A matrix, one should use rspde.make.A() in place of inla.spde.make.A(), and so on.
The main differences when comparing the arguments between the rSPDE-INLA implementation and the standard SPDE implementation in R-INLA, are the nu and rspde_order arguments, which are present in rSPDE-INLA implementation. We will see below how use these arguments.
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")
The observation matrix (A)
We now create the \(A\) matrix, that connects the mesh to the observation locations and then create the rSPDE model.
For this task, as we mentioned earlier, we need to use an rSPDEspecific function, whose name is very reminiscent to R-INLA’s standard SPDE approach, namely rspde.make.A() (in place of R-INLA’s inla.spde.make.A()). The reason for the need of this specific function is that the size of the \(A\) matrix depends on the order of the rational approximation. The details can be found in the introduction of the Rational approximation with the rSPDE package vignette.
The default order is 2 for our covariance-based rational approximation. As mentioned in the introduction of the Rational approximation with the rSPDE package vignette, an approximation of order 2 in the covariance-based rational approximation has approximately the same computational cost as the standard rational approximation of order 1.
Recall that the latent process \(u\) is a solution of \[(\kappa^2 I-\Delta)^{\alpha/2}(\tau u) = \mathcal{W},\] where \(\alpha = \nu + d/2\).
We want to estimate \(\tau,\kappa\) and \(\nu\).
There is also an option to fix the smoothness parameter \(\nu\) to some predefined value and only estimate \(\tau\) and \(\kappa\).
By default the rSPDE-INLA implementation estimates all the three parameters.
In this first example we will assume we want a rational approximation of order 2, where we also want to estimate the smoothness parameter \(\nu\).
To this end we can use the rspde.make.A() function. Since we will assume order 2 and that we want to estimate smoothness, which are the default options in this function, the required parameters are simply the mesh and the locations:
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:
Notice that this function is very reminiscent of R-INLA’s inla.spde2.matern() function.
This is a pattern we tried to keep consistent. All the rSPDE versions of some R-INLA function will either replace inla or inla.spde or inla.spde2 by rspde.
The inla.stack
Since the covariates are already evaluated at the observation locations, we only want to apply the \(A\) matrix to the spatial effect and not the fixed effects. We can use the inla.stack() function.
The difference, however, is that we need to use the function rspde.make.index() (in place of the standard inla.spde.make.index()) to create the index.
If one is using the default options, that is, to estimate the smoothness parameter \(\nu\) and to do a rational approximation of order 2, the usage of rspde.make.index() is identical to the usage of inla.spde.make.index():
We can then create the stack in a standard manner:
stk.dat <- inla.stack(
data = list(y = Y), A = list(Abar, 1), tag = "est",
effects = list(c(mesh.index,
list(Intercept = 1)),
list(long = inla.group(coords[, 1]),
lat = inla.group(coords[,2]),
seaDist = inla.group(seaDist))))Here the observation matrix \(A\) is applied to the spatial effect and the intercept while an identity observation matrix, denoted by \(1\), is applied to the covariates. This means the covariates are unaffected by the observation matrix.
The observation matrices in \(A=list(Abar,1)\) are used to link the corresponding elements in the effects-list to the observations. Thus in our model the latent spatial field mesh.index and the intercept are linked to the log-expectation of the observations, i.e. \(\eta(s)\), through the \(A\)-matrix. The covariates, on the other hand, are linked directly to \(\eta(s)\). The stk.dat object defined above implies the following principal linkage between model components and observations \[\eta(s) \sim A x(s) + A \text{ Intercept} + \text{long} + \text{lat}+ \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.
Here -1 is added to remove R’s implicit intercept, which is replaced by the explicit +Intercept from when we created the stack.
To fit the model we proceed as in the standard SPDE approach and we simply call inla(). It is noteworthy that we can have a performance enhancement by using the argument inla.mode = "experimental". As argued in R-INLA’s man page, not all features are available in inla.mode="experimental", and we will need to use the "classic" mode (the default mode) instead of "experimental" whenever we are doing predictions.
INLA 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):
##
## 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, ", " blas.num.threads = blas.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 = 3.81, Running = 24.6, Post = 0.0734, Total = 28.5
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 0.648 0.02 0.609 0.648 0.687 0.648 0
##
## Random effects:
## Name Model
## seaDist RW1 model
## field RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision parameter for the Gamma observations 13.12 0.900 11.444
## Precision for seaDist 9032.84 6385.059 2408.642
## Theta1 for field -1.25 0.497 -2.234
## Theta2 for field 1.09 0.251 0.599
## Theta3 for field -0.82 0.342 -1.488
## 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 13.085 14.981 13.019
## Precision for seaDist 7298.453 25848.274 5034.889
## Theta1 for field -1.244 -0.273 -1.235
## Theta2 for field 1.088 1.589 1.081
## Theta3 for field -0.823 -0.138 -0.834
##
## Marginal log-Likelihood: -1262.43
## 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 will be discussed later in this 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 will also be discussed later in this vignette.
rSPDE-INLA results
We can obtain outputs with respect to parameters in the original scale by using the function rspde.result():
## mean sd 0.025quant 0.5quant 0.975quant mode
## tau 0.324522 0.168490 0.108411 0.288546 0.760979 0.225534
## kappa 3.064570 0.779383 1.829630 2.968680 4.894780 2.775470
## nu 1.239940 0.285905 0.741572 1.220750 1.860740 1.171030We can also plot the posterior densities:
This function is reminiscent to the inla.spde.result() function with the main difference that it has the summary() and plot() methods implemented.
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 rspde.mesh.projector() function. This function has the same arguments as the function inla.mesh.projector(), with 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:
nxy <- c(150, 100)
projgrid <- rspde.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.
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")
Now, there are a few ways we could calculate the kriging prediction. The simplest way is to evaluate the mean of all individual random effects in the linear predictor and then to calculate the exponential of their sum (since \(\mu(s)=\exp(\eta(s))\) ). A more accurate way is to calculate the prediction jointly with the estimation, which unfortunately is quite computationally expensive if we do prediction on a fine grid. However, in this illustration, we proceed with this option to show how one can do it.
To this end, first, link the prediction coordinates to the mesh nodes through an \(A\) matrix
Since we are using distance to the sea as a covariate, we also have to calculate this covariate for the prediction locations.
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, list(Intercept = 1)),
list(long = inla.group(coord.prd[,
1]), lat = inla.group(coord.prd[, 2]),
seaDist = inla.group(seaDist.prd)))
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)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. Since we are doing predictions, we should not use R-INLA’s experimental mode.
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"))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(stk.all, "prd")$data
sd.prd <- m.prd <- matrix(NA, nxy[1], nxy[2])
m.prd[xy.in] <- rspde_fitprd$summary.fitted.values$mean[id.prd]
sd.prd[xy.in] <- rspde_fitprd$summary.fitted.values$sd[id.prd]Finally, we plot the results:
#Plot the predictions
grid.arrange(levelplot(m.prd, col.regions = tim.colors(99),
xlab = "", ylab = "", main = "mean",
scales = list(draw = FALSE)),
levelplot(sd.prd, col.regions = topo.colors(99),
xlab = "", ylab = "", scales = list(draw = FALSE),
main = "standard deviation"))
