Introduction

This document was written to explain how to build univariate spatiotemporal model with discrete time using INLA.

Note that INLA can be used with a range of data type because many different likelihood and link functions have been implemented in INLA. The different likelihood functions already implemented in INLA can be found by typing INLA::inla.list.models("likelihood") while the link functions can be found by typing INLA::inla.list.models("link").

To illustrate how these models are constructed, sugar maple data from Eastern North America will be used.

Load R package for the analysis

library(INLA)
library(raster)
library(sf)

Download and load the data

sugarRaw <- readRDS("assets/data/tree.rds")
climate <- readRDS("assets/data/climate_Present.rds")
region <- readRDS("assets/data/region.rds")

Organize environmental data

# Focus only on the data between 2010 and 2012
sugar <- sugarRaw[which(sugarRaw$yr > 2007),]

# Extract climate value for samples
climateSugar <- extract(climate, st_coordinates(sugar))
climate2007 <- scale(climateSugar[,1:2])

# Number of years
nyear <- length(unique(sugar$yr))

Building a rough polygon of the study area

# Extract extent of the sf data
extentRegion <- extent(sugar)

# Build polygon
polyRegion <- as(extentRegion, 'SpatialPolygons')

Build the INLA spatio-temporal model

Step 1 - Building the mesh

The first step that needs to be carried out to construct a spatiotemporal model is to build a Delaunay triangulation, a mesh, that will serve as a basis to construct the spatial component of the model through an SPDE approach. Properly choosing the structure of the mesh is essential for the model estimation to be carried out properly. Detailed explanation on the dos and don’ts of constructing a mesh is presented in section 1.3 of the R-INLA tutorial on SPDE models. For a quick overview, take a look at the mesh practical. For our illustration, let’s use the following mesh.

Spatial mesh

MeshSpace <- inla.mesh.2d(boundary = region, max.edge = 1,
                          min.angle = 20,
                          cutoff = 0.5,
                          offset = c(0.5,0.5),
                          crs = crs(region))

par(mar = c(1,1,1,1))

plot(MeshSpace, main = "", asp = 1)

# Number of knots in the mesh
MeshSpace$n

## [1] 523

Step 2 - Define the stochastic partial differential equation object

This is important to define the structure of the field on which the spatial contrain is imposed. In some ways, this serves as our “model variogram”. Recall that in INLA, it is the Matérn model variogram that is used. In this respect, in addition of the mesh, we have to define the scaling parameter (alpha), the prior range (prior.range) and the prior variance (prior.sigma).

alpha: This value needs to range between 0 and 2
prior.range: This is a vector of length 2 specifying the range (first) and the uncertainty around this range (second). Note that the uncertainty around the range needs to be smaller that the range itself.
prior.sigma: This is a vector of length 2 specifying the sill (first) and the uncertainty around this sill (second). Note that the uncertainty around the sill needs to be smaller that the range itself.

SPDE <- inla.spde2.pcmatern(mesh=MeshSpace,
                            alpha=2,
                            prior.range=c(0.1, 0.01),
                            prior.sigma=c(1, 0.01))

Step 3 - Priors and hyperparameters for temporal autocorrelation

hSpec <- list(theta=list(prior='pccor1', param=c(0, 0.9)))
precPrior <- list(prior='pc.prec', param=c(1, 0.01))

Step 4 - Index matrix

Here, we construct an “index matrix” used to extract the values of the latent spatial field at the observation locations.

Field <- inla.spde.make.index("field", n.spde=SPDE$n.spde,
                              n.group = nyear)

Step 5 - \(A\) matrix

We need to construct an \(A\) matrix. Recall from the lecture that the \(A\) matrix is sometimes referred to as the observation matrix and includes the parameters used to construct model. The discrete temporal groups are defined here.

# For estimation
A<-inla.spde.make.A(MeshSpace,
                    loc=st_coordinates(sugar),
                    group = sugar$yr - 2007)

Important note : The groups need to start at 1.

Step 6 - Organise the \(A\) matrix into a list

Alist <- as.list(rep(23,4))
Alist[[1]] <- A

Step 7 - Organise the effects (spatial autocorrelation structure and explanatory variables)

effect<- list(Field,
              bio1 =  climate2007[,1],
              bio2 = climate2007[,2],
              Intercept = rep(1,nrow(climate2007)))

Step 8 - Build `stack`

Here, we build a stack that will be passed to the inla function. Note that this is a special version of stack especially designed for the INLA package.

Stack <- inla.stack(data=list(tree = sugar$n_occ),
                    A = Alist,
                    effects = effect,
                    tag="basis")

Step 9 - Building the model (Finally !)

form <- tree ~ 0 + bio1 + bio2 + Intercept +
        f(field, # Space
          model=SPDE, # Space
          group=field.group, # Time 
          control.group=list(model='ar1', hyper=hSpec)) # Time 

model <- inla(form,
              data = inla.stack.data(Stack),
              family="nbinomial",
              control.family =list(link="log", hyper=list(theta=precPrior)),
              control.predictor=list(A=inla.stack.A(Stack),
                                                         compute=TRUE, link = 1),
              control.compute=list(waic=TRUE))

Important note : the argument group of the function f() refers to group defined in the object Field. As such, the “name” given in this object needs to be the name of the “field”.

Study the model output

summary(model)

## 
## Call:
##    c("inla(formula = form, family = \"nbinomial\", data = 
##    inla.stack.data(Stack), ", " control.compute = list(waic = TRUE), 
##    control.predictor = list(A = inla.stack.A(Stack), ", " compute = TRUE, 
##    link = 1), control.family = list(link = \"log\", ", " hyper = 
##    list(theta = precPrior)))") 
## Time used:
##     Pre = 1.92, Running = 183, Post = 0.815, Total = 185 
## Fixed effects:
##             mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## bio1      -0.009 0.001     -0.012   -0.009     -0.007 -0.009   0
## bio2      -0.005 0.001     -0.008   -0.005     -0.002 -0.005   0
## Intercept  0.087 0.005      0.077    0.087      0.096  0.087   0
## 
## Random effects:
##   Name     Model
##     field SPDE2 model
## 
## Model hyperparameters:
##                                                         mean    sd 0.025quant
## size for the nbinomial observations (1/overdispersion) 1.748 0.053      1.647
## Range for field                                        2.355 0.921      1.079
## Stdev for field                                        0.248 0.041      0.177
## GroupRho for field                                     0.971 0.030      0.889
##                                                        0.5quant 0.975quant
## size for the nbinomial observations (1/overdispersion)    1.747      1.856
## Range for field                                           2.180      4.642
## Stdev for field                                           0.245      0.337
## GroupRho for field                                        0.981      0.999
##                                                         mode
## size for the nbinomial observations (1/overdispersion) 1.744
## Range for field                                        1.878
## Stdev for field                                        0.238
## GroupRho for field                                     0.997
## 
## Expected number of effective parameters(stdev): 76.12(11.23)
## Number of equivalent replicates : 43.68 
## 
## Watanabe-Akaike information criterion (WAIC) ...: 20578.59
## Effective number of parameters .................: 79.96
## 
## Marginal log-Likelihood:  -10339.55 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

Plot prediction maps

# Dimention of the raster
rasterDim <- dim(climate)[1:2]

# Define basis of the map
mapBasis <- inla.mesh.projector(MeshSpace,
                               dims = rasterDim,
                               xlim = c(xmin(region), xmax(region)),
                               ylim = c(ymin(region), ymax(region)),
                               crs = crs(region))


# Calculate prediction
mapMean <- vector("list", length = nyear)
map.025 <- vector("list", length = nyear)
map.975 <- vector("list", length = nyear)

for(i in 1:nyear){
  mapMean[[i]] <- inla.mesh.project(mapBasis, 
                               model$summary.random$field$mean[Field$field.group == i])
  map.025[[i]] <- inla.mesh.project(mapBasis, 
                               model$summary.random$field$`0.025quant`[Field$field.group == i])
  map.975[[i]] <- inla.mesh.project(mapBasis, 
                               model$summary.random$field$`0.975quant`[Field$field.group == i])
}
  
# Transform map into a raster
rasterMean <- brick(climate, nl = 5)
raster.025 <- brick(climate, nl = 5)
raster.975 <- brick(climate, nl = 5)

for(i in 1:nyear){
  values(rasterMean)[,i] <- exp(as.vector(t(mapMean[[i]])))
  values(raster.025)[,i] <- exp(as.vector(t(map.025[[i]])))
  values(raster.975)[,i] <- exp(as.vector(t(map.975[[i]])))
}

# Mask the region
rasterMeanMask <- mask(rasterMean, region)
raster.025Mask <- mask(raster.025, region)
raster.975Mask <- mask(raster.975, region)

# Plot the results
for(i in 1:nyear){
par(mfrow = c(1,3))
plot(raster.025Mask[[i]], zlim = range(values(raster.025Mask[[i]]),
                              values(rasterMeanMask[[i]]),
                              values(raster.975Mask[[i]]), na.rm = TRUE),
     main = "2.5%")

plot(rasterMeanMask[[i]], zlim = range(values(raster.025Mask[[i]]),
                              values(rasterMeanMask[[i]]),
                              values(raster.975Mask[[i]]), na.rm = TRUE),
     main = "Mean")

plot(raster.975Mask[[i]], zlim = range(values(raster.025Mask[[i]]),
                              values(rasterMeanMask[[i]]),
                              values(raster.975Mask[[i]]), na.rm = TRUE),
     main = "97.5%")
}

Spatiotemporal model with INLA assuming discrete time

Guillaume Blanchet & Steve Vissault

November 21th, 2019

Introduction

Load R package for the analysis

Download and load the data

Organize environmental data

Building a rough polygon of the study area

Build the INLA spatio-temporal model

Step 1 - Building the mesh

Spatial mesh

Step 2 - Define the stochastic partial differential equation object

Step 3 - Priors and hyperparameters for temporal autocorrelation

Step 4 - Index matrix

Step 5 - \(A\) matrix

Step 6 - Organise the \(A\) matrix into a list

Step 7 - Organise the effects (spatial autocorrelation structure and explanatory variables)

Step 8 - Build `stack`

Step 9 - Building the model (Finally !)

Study the model output

Plot prediction maps

Spatiotemporal model with INLA assuming discrete time

Guillaume Blanchet & Steve Vissault

November 21th, 2019

Introduction

Load R package for the analysis

Download and load the data

Organize environmental data

Building a rough polygon of the study area

Build the INLA spatio-temporal model

Step 1 - Building the mesh

Spatial mesh

Step 2 - Define the stochastic partial differential equation object

Step 3 - Priors and hyperparameters for temporal autocorrelation

Step 4 - Index matrix

Step 5 - \(A\) matrix

Step 6 - Organise the \(A\) matrix into a list

Step 7 - Organise the effects (spatial autocorrelation structure and explanatory variables)

Step 8 - Build stack

Step 9 - Building the model (Finally !)

Study the model output

Plot prediction maps

Step 8 - Build `stack`