Tutorial 12.13 - Spatial and spatio-temporal models with INLA

# Tutorial 12.13 - Spatial and spatio-temporal models with INLA

23 Oct 2017

## Overview

Introduce what spatial and spatio-temporal models are.. Include point reference data...

## Simulating spatial data

We can simulate spatial data if we assume that the data can be drawn from a spatial field in which observations have some sort of spatial dependency. Observations closer together in space are likely to be more similar to one another. That is the similarity between locations is a function of the distance between the locations.

Lets start by generating 200 locations (quadrats, sites, etc) within a square domain. We will generate the coordinates of the locations such that the points are randomly scattered throughout the domain, however, we could alternatively sample in a way that yields a higher density of points in a particular region. Furthermore, rather than assume a perfectly square domain, we could sample from a polygon...

## generate points from a simple square domain
set.seed(1)
n <- 100
coords <- cbind(long=sample(1:n), lat=sample(1:n))
plot(coords)
coords.df <- data.frame(coords)
library(ggplot2)
ggplot(coords.df, aes(y=lat, x=long)) + geom_point()
set.seed(123)
n <- 100
coords1 <- cbind(long = sample(1:n/n - 0.5/n)^2, lat = sample(1:n/n - 0.5/n)^2)
plot(coords1, asp = 1)
set.seed(1)
library(sp)
## Create a spatialpolygons object - lets make use of the meuse.area matrix in the
## sp package
data(meuse.area)
meuse.sp <- SpatialPolygons(list(Polygons(list(Polygon(meuse.area)), "x")))
plot(meuse.sp)
coords.sp <- spsample(meuse.sp, n = 100, type = "random")
plot(coords.sp)

The spatial dependency between points ($s_i$ and $s_j$) could be defined by a Matern covariance structure. A Matern correlation function is: $$Cor_M(X(s_i), X(s_j)) = \frac{2^{1-\nu}}{\Gamma(\nu)}(\kappa\parallel s_i - s_j\parallel)^\nu K_\nu(\kappa\parallel s_i - s_j\parallel)$$ where $\parallel s_i - s_j\parallel$ is the euclidean distance between two points, $K_\nu$ is the modified Bessel function of the second order, $\kappa$ (kappa) is the scale parameter and $\nu$ (nu) is the smoothness parameter. The Matern covariance function is thus: $$\Sigma_{i,j}=\sigma_x Cor_M(X(s_i), X(s_j))$$ where $\sigma_x$ is the marginal variance of the process. This defines the process of contamination between points over the spatial domain. Therefore, the relative value of any point within this domain can be generated from a multivariate normal distribution with a mean of 0 and a variance of $\Sigma$.

kappa <- 7
nu <- 1
dmat <- dist(coords)
cor_m <- as.matrix((2^(1 - nu)/gamma(nu)) * ((kappa * dmat)^nu) * (besselK(dmat *
kappa, nu)))
cor_m[1:5, 1:5]
1             2             3             4             5
1  0.000000e+00  1.794100e-98 6.094506e-149 5.129406e-210  6.902439e-26
2  1.794100e-98  0.000000e+00  5.051134e-65 2.978713e-245  6.755479e-94
3 6.094506e-149  5.051134e-65  0.000000e+00 3.006513e-233 4.077963e-152
4 5.129406e-210 2.978713e-245 3.006513e-233  0.000000e+00 7.386427e-236
5  6.902439e-26  6.755479e-94 4.077963e-152 7.386427e-236  0.000000e+00

In generating observed responses at the locations, it becomes realistic to add additional stochasticity due to measurement error. So that the value of the observed response at locations $s_i$ is equal to the realization of the underlying spatial field (process) at that location ($x(s_i)$) plus some additional measurement error ($e_i$) that we assume to be independent and normally distributed: $$y(s_i) = x(s_i) + e_i \hspace{1cm} e_i \sim{}N(0, \sigma^2_e)$$ where $\sigma^2_e$ is the 'nugget' effect. Thus the marginal distribution of the observed data ($y(s)$) is: $\sigma^2_eI + \Sigma$

beta0 <- 10
sigma2e <- 0.3
sigma2x <- 5

diag(cor_m) <- 1
cov_m <- sigma2e * diag(n) + sigma2x * cor_m

L <- chol(cov_m)
set.seed(234)
y <- beta0 + drop(rnorm(n) %*% L)

coords.df <- data.frame(coords.df, y = y)
long lat
[1,]   27  66
[2,]   37  35
[3,]   57  27
[4,]   89  97
[5,]   20  61
[6,]   86  21
library(ggplot2)
ggplot(coords.df, aes(y = lat, x = long)) + geom_point(size = y/2)

## SPDE

Give an intro to stochastic partial differential equations (SPDE) and the computational advantages.

### Mesh

#### Convex hull based meshes

The inla.mesh.2d() function can be used to generate a Constrained Refined Delaunay Triangulation (CRDT). To demonstrate the various mesh generation options, we will generate a number of meshes. As there are potentially more options relevant to spatial points within an irregular polygon domain rather than a perfectly square domain, I will first demonstrate meshes for the coords.sp coordinates based on points sampled from the meuse.area object earlier.

The function must be supplied with either the coordinates of locations or else the location domain (coordinates of a hull around all the points) in which to generate the triangles. In addition, we must indicate the minimum length of triangle edges. In doing so, we indirectly dictate the granularity with which the domain is broken up into triangles (the smaller the edges, the smaller the triangles, and thus the more of them).

As the variance properties of triangles around the edge are half that of triangles within the interior, it is a good idea to define an outer edge as a buffer.

ArgumentDescription
loc a matrix of coordinates used as initial triangulation nodes
loc.domain coordinates of a polygon that defines the spatial domain
max.edge maximum permissible triangle edge length for the inner (and outer) region. The value(s) should be on a scale relevant to the magnitude of the coordinates. The lower the values, the more triangles.
offset an amount to expand the distance between the points and the edge of the inner (and outer) region. Positive values are treated as absolute distances and negative numbers are treated as multipliers.
cuttoff the minimum distance allowed between points. This allows points that are very close together to be placed in the same triangle. Of particular concern, is that we don't want triangles with very sharp angles as these will be less effective when it comes to projections.
library(INLA)
## Convert spatialpolygon object into matrix
coords.sp <- as.matrix(as.data.frame(coords.sp))
mesh1 <- inla.mesh.2d(loc = coords.sp, max.edge = c(2000, 2000))
plot(mesh1)
points(coords.sp, col = "red", pch = 16)
mesh2 <- inla.mesh.2d(loc = coords.sp, max.edge = c(500, 2000))
plot(mesh2)
points(coords.sp, col = "red", pch = 16)
mesh3 <- inla.mesh.2d(loc = coords.sp, max.edge = c(500, 2000), offset = c(1000,
3000))
plot(mesh3)
points(coords.sp, col = "red", pch = 16)
mesh4 <- inla.mesh.2d(loc = coords.sp, max.edge = c(500, 2000), offset = c(-0.2,
-0.5))
plot(mesh4)
points(coords.sp, col = "red", pch = 16)
mesh5 <- inla.mesh.2d(loc = coords.sp, max.edge = c(2000, 2000), cutoff = 500)
plot(mesh5)
points(coords.sp, col = "red", pch = 16)
mesh6 <- inla.mesh.2d(loc = coords.sp, max.edge = c(500, 2000), boundary = list(inla.mesh.segment(meuse.area)))
plot(mesh6)
points(coords.sp, col = "red", pch = 16)
lines(meuse.area, lwd = 2)
## or if we just have a matrix of coordinates around the spatial domain
mesh7 <- inla.mesh.2d(loc.domain = meuse.area, max.edge = c(500, 2000))
plot(mesh7)

#### Non-convex hull meshes

Alternatively, we can base the domain on a non-convex hull by using the inla.nonconvex.hull() function.

library(INLA)
## Convert spatialpolygon object into matrix
coords.sp <- as.matrix(as.data.frame(coords.sp))
boundary <- inla.nonconvex.hull(points = coords.sp)
mesh8 <- inla.mesh.2d(boundary = boundary, max.edge = c(500, 2000))
plot(mesh8)
points(coords.sp, col = "red", pch = 16)

Returning to our square domain, we want to establish a mesh that has a decent buffer (outer edge)

mesh <- inla.mesh.2d(coords, max.edge = c(20, 40))
plot(mesh)
points(coords, col = "red", pch = 16)

This mesh is a list containing:

str(mesh)
List of 8
$meta :List of 5 ..$ call        : language inla.mesh.create(loc = rbind(loc, mesh2$loc), boundary = boundary[[2]], interior = c(boundary2, interior2), extend = list(n = n[2], offset = offset[2]), ... ..$ fmesher.args: chr "--input=input.s --cutoff=1e-12  --interior=input.segm.int.idx  --interiorgrp=input.segm.int.grp --cet=16,-0.15 --rcdt=21,40,40 "| __truncated__
..$time : num [1:5, 1:5] 0 0 0 0 0 0 0 0 0 0 ... .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : chr [1:5] "pre" "fmesher" "post" "object" ...
.. .. ..$: NULL ..$ prefix      : chr "/tmp/RtmpnLy1Nt/fmesher3e541a5fea57."
..$is.refined : logi TRUE$ manifold: chr "R2"
$n : int 223$ loc     : num [1:223, 1:3] 8.28 90.72 106.58 106.58 100.72 ...
$graph :List of 5 ..$ tv : int [1:415, 1:3] 29 127 73 64 200 26 50 49 12 23 ...
..$vt : int [1:223, 1] 12 21 36 414 53 7 10 85 11 35 ... ..$ tt : int [1:415, 1:3] 26 212 161 91 53 372 69 113 65 111 ...
..$tti: int [1:415, 1:3] 1 3 1 3 1 1 2 1 2 2 ... ..$ vv :Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
.. .. ..@ i       : int [1:1274] 26 37 61 95 192 207 20 30 55 127 ...
.. .. ..@ j       : int [1:1274] 0 0 0 0 0 0 1 1 1 1 ...
.. .. ..@ Dim     : int [1:2] 223 223
.. .. ..@ Dimnames:List of 2
.. .. .. ..$: NULL .. .. .. ..$ : NULL
.. .. ..@ x       : num [1:1274] 1 1 1 1 1 1 1 1 1 1 ...
.. .. ..@ factors : list()
$segm :List of 2 ..$ bnd:List of 5
.. ..$loc : NULL .. ..$ idx   : int [1:29, 1:2] 193 218 211 217 194 195 196 197 221 213 ...
.. ..$grp : int [1:29, 1] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ is.bnd: logi TRUE
.. ..$crs : NULL .. ..- attr(*, "class")= chr "inla.mesh.segment" ..$ int:List of 5
.. ..$loc : NULL .. ..$ idx   : int [1:38, 1:2] 1 38 17 25 9 33 14 31 2 21 ...
.. ..$grp : int [1:38, 1] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ is.bnd: logi FALSE
.. ..$crs : NULL .. ..- attr(*, "class")= chr "inla.mesh.segment"$ idx     :List of 2
..$loc : int [1:100] 39 40 41 42 43 44 45 46 47 48 ... ..$ lattice: NULL
$crs : NULL - attr(*, "class")= chr "inla.mesh" ## The number of vertices in the mesh mesh$n
[1] 223
## The indices of the vertices that correspond to the coordinates of the observed
## locations
mesh$idx$loc
[1]  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56
[19]  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74
[37]  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92
[55]  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108 109 110
[73] 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
[91] 129 130 131 132 133 134 135 136 137 138

### SPDE model

spde <- inla.spde2.matern(mesh, alpha = 2)
str(spde)
List of 7
$model : chr "matern"$ n.spde    : int 223
$n.theta : int 2$ param.inla:List of 15
..$n : int 223 ..$ n.theta      : int 2
..$M0 :Formal class 'dgTMatrix' [package "Matrix"] with 6 slots .. .. ..@ i : int [1:223] 0 1 2 3 4 5 6 7 8 9 ... .. .. ..@ j : int [1:223] 0 1 2 3 4 5 6 7 8 9 ... .. .. ..@ Dim : int [1:2] 223 223 .. .. ..@ Dimnames:List of 2 .. .. .. ..$ : NULL
.. .. .. ..$: NULL .. .. ..@ x : num [1:223] 190 171 208 117 129 ... .. .. ..@ factors : list() ..$ M1           :Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
.. .. ..@ i       : int [1:1497] 0 26 37 61 95 192 207 1 20 30 ...
.. .. ..@ j       : int [1:1497] 0 0 0 0 0 0 0 1 1 1 ...
.. .. ..@ Dim     : int [1:2] 223 223
.. .. ..@ Dimnames:List of 2
.. .. .. ..$: NULL .. .. .. ..$ : NULL
.. .. ..@ x       : num [1:1497] 5.1604 -2.1409 -1.9881 -0.0821 -0.0883 ...
.. .. ..@ factors : list()
..$M2 :Formal class 'dgTMatrix' [package "Matrix"] with 6 slots .. .. ..@ i : int [1:4041] 0 7 16 26 37 61 65 84 95 108 ... .. .. ..@ j : int [1:4041] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..@ Dim : int [1:2] 223 223 .. .. ..@ Dimnames:List of 2 .. .. .. ..$ : NULL
.. .. .. ..$: NULL .. .. ..@ x : num [1:4041] 0.1841 0.0137 0.0143 -0.1095 -0.1042 ... .. .. ..@ factors : list() ..$ B0           : num [1:223, 1:3] 0 0 0 0 0 0 0 0 0 0 ...
..$B1 : num [1:223, 1:3] 0 0 0 0 0 0 0 0 0 0 ... ..$ B2           : num [1:223, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$: NULL .. .. ..$ : chr [1:3] "B" "" ""
..$BLC : num [1:6, 1:3] 0 0 0 0 -2.53 ... .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : chr [1:6] "theta.1" "theta.2" "tau.1" "kappa.1" ...
.. .. ..$: NULL ..$ theta.mu     : num [1:2] 1.13 -2.4
..$theta.Q : num [1:2, 1:2] 0.1 0 0 0.1 ..$ transform    : chr "identity"
..$theta.initial: num [1:2] 1.13 -2.4 ..$ fixed        : logi [1:2] FALSE FALSE
..$theta.fixed : num(0)$ f         :List of 4
..$model : chr "spde2" ..$ n              : int 223
..$spde2.transform: chr "identity" ..$ hyper.default  :List of 2
.. ..$theta1:List of 3 .. .. ..$ prior  : chr "mvnorm"
.. .. ..$param : num [1:6] 1.13 -2.4 0.1 0 0 ... .. .. ..$ initial: num 1.13
.. ..$theta2:List of 1 .. .. ..$ initial: num -2.4
$BLC : num [1:6, 1:3] 0 0 0 0 -2.53 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:6] "theta.1" "theta.2" "tau.1" "kappa.1" ...
.. ..$: NULL$ mesh      :List of 8
..$meta :List of 5 .. ..$ call        : language inla.mesh.create(loc = rbind(loc, mesh2$loc), boundary = boundary[[2]], interior = c(boundary2, interior2), extend = list(n = n[2], offset = offset[2]), ... .. ..$ fmesher.args: chr "--input=input.s --cutoff=1e-12  --interior=input.segm.int.idx  --interiorgrp=input.segm.int.grp --cet=16,-0.15 --rcdt=21,40,40 "| __truncated__
.. ..$time : num [1:5, 1:5] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..- attr(*, "dimnames")=List of 2 .. .. .. ..$ : chr [1:5] "pre" "fmesher" "post" "object" ...
.. .. .. ..$: NULL .. ..$ prefix      : chr "/tmp/RtmpnLy1Nt/fmesher3e541a5fea57."
.. ..$is.refined : logi TRUE ..$ manifold: chr "R2"
..$n : int 223 ..$ loc     : num [1:223, 1:3] 8.28 90.72 106.58 106.58 100.72 ...
..$graph :List of 5 .. ..$ tv : int [1:415, 1:3] 29 127 73 64 200 26 50 49 12 23 ...
.. ..$vt : int [1:223, 1] 12 21 36 414 53 7 10 85 11 35 ... .. ..$ tt : int [1:415, 1:3] 26 212 161 91 53 372 69 113 65 111 ...
.. ..$tti: int [1:415, 1:3] 1 3 1 3 1 1 2 1 2 2 ... .. ..$ vv :Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
.. .. .. ..@ i       : int [1:1274] 26 37 61 95 192 207 20 30 55 127 ...
.. .. .. ..@ j       : int [1:1274] 0 0 0 0 0 0 1 1 1 1 ...
.. .. .. ..@ Dim     : int [1:2] 223 223
.. .. .. ..@ Dimnames:List of 2
.. .. .. .. ..$: NULL .. .. .. .. ..$ : NULL
.. .. .. ..@ x       : num [1:1274] 1 1 1 1 1 1 1 1 1 1 ...
.. .. .. ..@ factors : list()
..$segm :List of 2 .. ..$ bnd:List of 5
.. .. ..$loc : NULL .. .. ..$ idx   : int [1:29, 1:2] 193 218 211 217 194 195 196 197 221 213 ...
.. .. ..$grp : int [1:29, 1] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..$ is.bnd: logi TRUE
.. .. ..$crs : NULL .. .. ..- attr(*, "class")= chr "inla.mesh.segment" .. ..$ int:List of 5
.. .. ..$loc : NULL .. .. ..$ idx   : int [1:38, 1:2] 1 38 17 25 9 33 14 31 2 21 ...
.. .. ..$grp : int [1:38, 1] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..$ is.bnd: logi FALSE
.. .. ..$crs : NULL .. .. ..- attr(*, "class")= chr "inla.mesh.segment" ..$ idx     :List of 2
.. ..$loc : int [1:100] 39 40 41 42 43 44 45 46 47 48 ... .. ..$ lattice: NULL
..$crs : NULL ..- attr(*, "class")= chr "inla.mesh" - attr(*, "class")= chr [1:3] "inla.spde2" "inla.spde" "inla.model.class" The SPDE model can get quite complex. So to assist with keeping track of which elements relate to what effects, we can create an index. s.index <- inla.spde.make.index("spatial.field", n.spde = spde$n.spde)
str(s.index)
List of 3
$spatial.field : int [1:223] 1 2 3 4 5 6 7 8 9 10 ...$ spatial.field.group: int [1:223] 1 1 1 1 1 1 1 1 1 1 ...
$spatial.field.repl : int [1:223] 1 1 1 1 1 1 1 1 1 1 ... The vector called spatial.field, contains the indices for the spatial field. In this case, our spatial data are all in the one group. ## The projector matrix The SPDE model is defined on a mesh with$m$dimensions, whereas the response is defined at$n$locations. We need a way to link the$m$mesh vertices to the$n$responses. This is achieved via a projector matrix ($\mathbf{A}$). This projector matrix is build using the inla.spde.make.A() function. The value of any point is created by projecting this point onto a plane formed by the triangle it is positioned in. The exact value of the point is the weighted average of three weights associated with each vertex. To illustrate this, it might be best to deviate from our working set of coordinates and mesh and generate a very small number of points. To better appreciate the weights that are applied to each vertex, we will create two examples. One where the points are on the vertices (and thus weights are either 0 or one) and another example where the points are not necessarily on the vertices and thus each point will be associated with up to three non-zero weights (that must collectively add up to 1). set.seed(12) n <- 4 coords2 <- cbind(long = sample(1:n), lat = sample(1:n)) plot(coords2) mesh9 <- inla.mesh.2d(coords2, max.edge = c(50, 50)) plot(mesh9) points(coords2, col = "red", pch = 16) A9 <- inla.spde.make.A(mesh9, loc = coords2) str(A9) Formal class 'dgCMatrix' [package "Matrix"] with 6 slots ..@ i : int [1:12] 2 1 3 0 2 0 1 2 3 1 ... ..@ p : int [1:63] 0 0 0 0 0 0 0 0 0 1 ... ..@ Dim : int [1:2] 4 62 ..@ Dimnames:List of 2 .. ..$ : NULL
.. ..$: NULL ..@ x : num [1:12] 0 0 0 0 0 1 1 1 1 0 ... ..@ factors : list() dim(A9) [1] 4 62 A9 4 x 62 sparse Matrix of class "dgCMatrix" [1,] . . . . . . . . . . . . . . . . . . . . 0 . . . . 1 . . . . . . . . 0 . . [2,] . . . . . . . . . . . . . . . . 0 . . . . . . . . . 1 . . . . . 0 . . . . [3,] . . . . . . . . 0 . . . . . . . . . . . . . . 0 . . . 1 . . . . . . . . . [4,] . . . . . . . . . . . . . . . . . . . 0 . . . . . . . . 1 . . . . 0 . . . [1,] . . . . . . . . . . . . . . . . . . . . . . . . . [2,] . . . . . . . . . . . . . . . . . . . . . . . . . [3,] . . . . . . . . . . . . . . . . . . . . . . . . . [4,] . . . . . . . . . . . . . . . . . . . . . . . . . table(rowSums(A9 > 0)) 1 4 table(rowSums(A9)) 1 4 table(colSums(A9) > 0) FALSE TRUE 58 4 The A9 object is a sparse matrix that relates (via weights) the observed points (rows) to all the vertices (columns). The mesh contained 62 vertices. However, of these, only 4 columns contain data. set.seed(12) n <- 4 coords2 <- cbind(long = sample(1:n), lat = sample(1:n)) plot(coords2) boundary <- cbind(c(0, 5, 5, 0), c(0, 0, 5, 5)) mesh10 <- inla.mesh.2d(loc.domain = boundary, max.edge = c(4, 4)) plot(mesh10) points(coords2, col = "red", pch = 16) A10 <- inla.spde.make.A(mesh10, loc = coords2) str(A10) Formal class 'dgCMatrix' [package "Matrix"] with 6 slots ..@ i : int [1:12] 1 3 0 0 1 2 3 1 2 2 ... ..@ p : int [1:67] 0 0 0 0 0 0 0 0 0 1 ... ..@ Dim : int [1:2] 4 66 ..@ Dimnames:List of 2 .. ..$ : NULL
.. ..$: NULL ..@ x : num [1:12] 0.382 0.349 0.197 0.185 0.539 ... ..@ factors : list() dim(A10) [1] 4 66 round(A10, 2) 4 x 66 sparse Matrix of class "dgCMatrix" [1,] . . . . . . . . . . . . . 0.2 . . . . . . . . . . . . 0.18 . . [2,] . . . . . . . . 0.38 . . . . . . . . . . . . . . . . . . . 0.54 [3,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.20 [4,] . . . . . . . . . 0.35 . . . . . . . . . . . . . . . . . . 0.13 [1,] . . 0.62 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [2,] 0.08 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [3,] 0.71 0.09 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [4,] . 0.52 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [1,] . . . . [2,] . . . . [3,] . . . . [4,] . . . . table(rowSums(A10 > 0)) 3 4 table(rowSums(A10)) 1 4 table(colSums(A10) > 0) FALSE TRUE 58 8 We now return to our data.. A <- inla.spde.make.A(mesh, loc = coords) str(A) Formal class 'dgCMatrix' [package "Matrix"] with 6 slots ..@ i : int [1:300] 23 17 86 46 54 15 97 83 60 92 ... ..@ p : int [1:224] 0 1 2 2 3 3 3 3 5 7 ... ..@ Dim : int [1:2] 100 223 ..@ Dimnames:List of 2 .. ..$ : NULL
.. ..$: NULL ..@ x : num [1:300] 0 0 0 0 0 0 0 0 0 0 ... ..@ factors : list() dim(A) [1] 100 223 table(rowSums(A > 0)) 1 100 table(rowSums(A)) 1 100 table(colSums(A) > 0) FALSE TRUE 123 100 ## The stack It is possible for models to include multiple projector matrices - eg. Lindgren, 2012, Cameletti et al, 2012 anmd Lindgren and Rue, 2013. The inla stack allows us to work with predictors that include terms with different dimensions. stack <- inla.stack(data = list(y = coords.df$y), A = list(A, 1), effects = list(s.index,
list(Intercept = rep(1, nrow(coords.df)))), tag = "est")

Note that the inla.stack() function automatically eliminates all the elements associated with columns that sum to zero.

dim(inla.stack.A(stack))
[1] 100 101

## Fitting the INLA model

data.inla <- inla(y ~ 0 + Intercept + f(spatial.field, model = spde), data = inla.stack.data(stack),
control.predictor = list(A = inla.stack.A(stack)))
summary(data.inla)
Call:
c("inla(formula = y ~ 0 + Intercept + f(spatial.field, model = spde), ",  "    data = inla.stack.data(stack), control.predictor = list(A = inla.stack.A(stack)))" )

Time used:
Pre-processing    Running inla Post-processing           Total
0.5985          0.8630          0.0458          1.5073

Fixed effects:
mean     sd 0.025quant 0.5quant 0.975quant    mode kld
Intercept 10.008 0.2266     9.5619  10.0078     10.455 10.0073   0

Random effects:
Name	  Model
spatial.field   SPDE2 model

Model hyperparameters:
mean     sd 0.025quant 0.5quant
Precision for the Gaussian observations  0.2197 0.0278     0.1649   0.2203
Theta1 for spatial.field                 2.1413 2.5020    -2.4934   2.0218
Theta2 for spatial.field                -0.5379 1.6146    -3.2479  -0.7125
0.975quant   mode
Precision for the Gaussian observations     0.2736  0.224
Theta1 for spatial.field                    7.3201  1.591
Theta2 for spatial.field                    3.0353 -1.363

Expected number of effective parameters(std dev): 3.502(2.897)
Number of equivalent replicates : 28.55

Marginal log-Likelihood:  -234.61

Since we have simulated these data, we can see how well we were able to recover the defining parameters:

• Intercept = 10
• sigma2e=0.3
• sigma2x=5
• kappa=7
• nu=1

## Intercept
data.inla$summary.fix mean sd 0.025quant 0.5quant 0.975quant mode Intercept 10.00803 0.2266212 9.561947 10.00776 10.45499 10.00725 kld Intercept 7.60943e-13 ## sigma2e - note, this is on a precision scale and so needs to be converted to ## standard deviation scale post.se <- inla.tmarginal(function(x) sqrt(1/x), data.inla$marginals.hy[[1]])

inla.emarginal(function(x) x, post.se)
[1] 2.146443
inla.hpdmarginal(0.95, post.se)
low     high
level:0.95 1.895175 2.431498
data.inla.field <- inla.spde2.result(data.inla, "spatial.field", spde, do.transf = TRUE)
inla.emarginal(function(x) x, data.inla.field$marginals.kappa[[1]]) [1] 2.786834 inla.emarginal(function(x) x, data.inla.field$marginals.variance.nominal[[1]])
[1] 0.04743584

### Prediction

Prediction involves projecting the fitted model into the mesh at the specific spatial locations of interest (i.e. our spatial grid). Recall that the spatial model is included as a random effect and therefore all coeffients sum to zero - that is, they are deviations from the overal mean (Intercept). Hence to project onto the response scale, we must add the estimated intercept.
long lat         y
1   27  66 11.521206
2   37  35  5.273678
3   57  27  6.548568
4   89  97 13.387033
5   20  61 13.359189
6   86  21 10.322624
coords.grid <- as.matrix(expand.grid(long = seq(0, 100, len = 100), lat = seq(0,
100, len = 100)))
long lat
[1,] 0.000000   0
[2,] 1.010101   0
[3,] 2.020202   0
[4,] 3.030303   0
[5,] 4.040404   0
[6,] 5.050505   0
data.inla.projector <- inla.mesh.projector(mesh, loc = coords.grid)
newdata <- data.frame(coords.grid, mean = inla.mesh.project(data.inla.projector,
data.inla$summary.random$spatial.field$mean) + data.inla$summary.fixed$mean[1]) str(newdata$mean)
num [1:10000] 9.98 9.98 9.98 9.98 9.97 ...
ggplot(newdata, aes(y = long, x = lat)) + geom_tile(aes(fill = mean))
ggplot(coords.df, aes(y = long, x = lat)) + geom_point(aes(color = y), size = 2)
library(mgcv)
dat.gamm <- gam(y ~ s(long, lat), data = coords.df)
newdata$y1 <- predict(dat.gamm, newdata = newdata) ggplot(newdata, aes(y = long, x = lat)) + geom_tile(aes(fill = y1)) coords.df$pred = data.inla$summary.fix[1, 1] + drop(A %*% data.inla$summary.random$spatial.field$mean)
ggplot(coords.df, aes(y = pred, x = y)) + geom_point()
cor(coords.df$y, coords.df$pred)
[1] 0.8134691
ggplot(coords.df, aes(y = predict(dat.gamm), x = y)) + geom_point()
cor(coords.df$y, predict(dat.gamm)) [1] 0.2471316 Now lets do the same, yet also include latitude and longitude as 'fixed' effects. To do so, we generate two stacks, one that represents the observed data, and another that represents the prediction data (with NA for the response). These stacks are then stacked together into a single stack. stack.e <- inla.stack(data = list(y = coords.df$y), A = list(A, 1, 1, 1), effects = list(s.index,
list(Intercept = rep(1, nrow(coords.df))), list(long = coords.df$long), list(lat = coords.df$lat)),
tag = "est")

newdata = data.frame(long = mesh$loc[, 1], lat = mesh$loc[, 2])

A.p = Diagonal(n = nrow(newdata))
stack.p = inla.stack(data = list(y = NA), A = list(A.p, 1, 1, 1), effects = list(s.index,
list(Intercept = rep(1, nrow(newdata))), list(lon = newdata$lon), list(lat = newdata$lat)),
tag = "pred")
stack = inla.stack(stack.e, stack.p)

Note that the inla.stack() function automatically eliminates all the elements associated with columns that sum to zero.

dim(inla.stack.A(stack))
[1] 323 716

The fitted model then includes fixed effects of lat and long. Note also, that since we are going to want to access the estimated fitted values of the prediction set, we must include the compute=TRUE argument to the control.predictor() parameter.
data.inla <- inla(y ~ 0 + Intercept + lat + long + f(spatial.field, model = spde),
data = inla.stack.data(stack), control.predictor = list(A = inla.stack.A(stack),
compute = TRUE))
summary(data.inla)
Call:
c("inla(formula = y ~ 0 + Intercept + lat + long + f(spatial.field, ",  "    model = spde), data = inla.stack.data(stack), control.predictor = list(A = inla.stack.A(stack), ",  "    compute = TRUE))")

Time used:
Pre-processing    Running inla Post-processing           Total
0.4937          1.5364          0.1319          2.1620

Fixed effects:
mean     sd 0.025quant 0.5quant 0.975quant   mode kld
Intercept 8.6925 0.5964     7.5101   8.6948     9.8605 8.6993   0
lat       0.0156 0.0076     0.0006   0.0156     0.0307 0.0156   0
long      0.0104 0.0076    -0.0046   0.0103     0.0255 0.0103   0

Random effects:
Name	  Model
spatial.field   SPDE2 model

Model hyperparameters:
mean     sd 0.025quant 0.5quant
Precision for the Gaussian observations  0.2312 0.0352     0.1609   0.2325
Theta1 for spatial.field                 3.4591 2.9245    -0.7670   2.9735
Theta2 for spatial.field                -0.0356 1.0841    -1.5648  -0.2240
0.975quant    mode
Precision for the Gaussian observations     0.2966  0.2392
Theta1 for spatial.field                   10.2423  0.8022
Theta2 for spatial.field                    2.4962 -1.0336

Expected number of effective parameters(std dev): 8.322(5.32)
Number of equivalent replicates : 12.02

Marginal log-Likelihood:  -250.03
Posterior marginals for linear predictor and fitted values computed
long lat         y      pred
1   27  66 11.521206 10.038941
2   37  35  5.273678  9.906943
3   57  27  6.548568  9.913135
4   89  97 13.387033 10.209700
5   20  61 13.359189 10.068761
6   86  21 10.322624  9.988433
coords.grid <- as.matrix(expand.grid(long = seq(0, 100, len = 100), lat = seq(0,
100, len = 100)))
long lat
[1,] 0.000000   0
[2,] 1.010101   0
[3,] 2.020202   0
[4,] 3.030303   0
[5,] 4.040404   0
[6,] 5.050505   0
data.inla.projector <- inla.mesh.projector(mesh, loc = coords.grid)

p.index = inla.stack.index(stack, "pred")$data data.means = data.inla$summary.fitted.values$mean[p.index] fit = sd = NULL newdata = cbind(coords.grid, mean = c(inla.mesh.project(data.inla.projector, data.means))) newdata = data.frame(newdata) ggplot(newdata, aes(y = long, x = lat)) + geom_tile(aes(fill = mean)) ggplot(coords.df, aes(y = long, x = lat)) + geom_point(aes(color = y), size = 2) library(mgcv) dat.gamm <- gam(y ~ s(long, lat), data = coords.df) newdata$y1 <- predict(dat.gamm, newdata = newdata)
ggplot(newdata, aes(y = long, x = lat)) + geom_tile(aes(fill = y1))
coords.df$pred = data.inla$summary.fix[1, 1] + drop(A %*% data.inla$summary.random$spatial.field$mean) ggplot(coords.df, aes(y = pred, x = y)) + geom_point() cor(coords.df$y, coords.df$pred) [1] 0.7852556 ggplot(coords.df, aes(y = predict(dat.gamm), x = y)) + geom_point() cor(coords.df$y, predict(dat.gamm))
[1] 0.2471316

## Spatio-temporal

Illustrate spatial pattern at discrete time intervals. We will start by just extending the previous spatial example. Lets say we had ten years of data.

### Simulate data

rspde <- function(coords, kappa, variance = 1, alpha = 2, n = 1, mesh, verbose = FALSE,
seed, return.attributes = FALSE) {
t0 <- Sys.time()
theta <- c(-0.5 * log(4 * pi * variance * kappa^2), log(kappa))
if (verbose)
cat("theta =", theta, "\n")
if (missing(mesh)) {
mesh.pars <- c(0.5, 1, 0.1, 0.5, 1) * sqrt(alpha - ncol(coords)/2)/kappa
if (verbose)
cat("mesh.pars =", mesh.pars, "\n")
attributes <- list(mesh = inla.mesh.2d(, coords[chull(coords), ], max.edge = mesh.pars[1:2],
cutoff = mesh.pars[3], offset = mesh.pars[4:5]))
if (verbose)
cat("n.mesh =", attributes$mesh$n, "\n")
} else attributes <- list(mesh = mesh)
attributes$spde <- inla.spde2.matern(attributes$mesh, alpha = alpha)
attributes$Q <- inla.spde2.precision(attributes$spde, theta = theta)
attributes$A <- inla.mesh.project(mesh = attributes$mesh, loc = coords)$A if (n == 1) result <- drop(attributes$A %*% inla.qsample(Q = attributes$Q, constr = attributes$spde$f$extraconstr))
t1 <- Sys.time()
result <- inla.qsample(n, attributes$Q, seed = ifelse(missing(seed), 0, seed), constr = attributes$spde$f$extraconstr)
if (nrow(result) < nrow(attributes$A)) { result <- rbind(result, matrix(NA, nrow(attributes$A) - nrow(result), ncol(result)))
dimnames(result)[[1]] <- paste("x", 1:nrow(result), sep = "")
for (j in 1:ncol(result)) result[, j] <- drop(attributes$A %*% result[1:ncol(attributes$A),
j])
} else {
for (j in 1:ncol(result)) result[1:nrow(attributes$A), j] <- drop(attributes$A %*%
result[, j])
result <- result[1:nrow(attributes$A), ] } t2 <- Sys.time() attributes$cpu <- c(prep = t1 - t0, sample = t2 - t1, total = t2 - t0)
if (return.attributes)
attributes(result) <- c(attributes(result), attributes)
return(drop(result))
}

We start by making $k=10$ realizations from the mesh we defined earlier.

k = 10
kappa = 7
nu = 1
rho = 0.7  # degree of temporal autocorrelation
x.k <- rspde(coords, kappa = 7, n = k, mesh = mesh, return.attributes = TRUE)
class(x.k)
[1] "matrix"
str(x.k)
num [1:100, 1:10] -0.000627 0.058355 0.020546 -0.183822 -0.055067 ...
- attr(*, "dimnames")=List of 2
..$: chr [1:100] "x1" "x2" "x3" "x4" ... ..$ : chr [1:10] "sample1" "sample2" "sample3" "sample4" ...
- attr(*, "mesh")=List of 8
..$meta :List of 5 .. ..$ call        : language inla.mesh.create(loc = rbind(loc, mesh2$loc), boundary = boundary[[2]], interior = c(boundary2, interior2), extend = list(n = n[2], offset = offset[2]), ... .. ..$ fmesher.args: chr "--input=input.s --cutoff=1e-12  --interior=input.segm.int.idx  --interiorgrp=input.segm.int.grp --cet=16,-0.15 --rcdt=21,40,40 "| __truncated__
.. ..$time : num [1:5, 1:5] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..- attr(*, "dimnames")=List of 2 .. .. .. ..$ : chr [1:5] "pre" "fmesher" "post" "object" ...
.. .. .. ..$: NULL .. ..$ prefix      : chr "/tmp/RtmpnLy1Nt/fmesher3e541a5fea57."
.. ..$is.refined : logi TRUE ..$ manifold: chr "R2"
..$n : int 223 ..$ loc     : num [1:223, 1:3] 8.28 90.72 106.58 106.58 100.72 ...
..$graph :List of 5 .. ..$ tv : int [1:415, 1:3] 29 127 73 64 200 26 50 49 12 23 ...
.. ..$vt : int [1:223, 1] 12 21 36 414 53 7 10 85 11 35 ... .. ..$ tt : int [1:415, 1:3] 26 212 161 91 53 372 69 113 65 111 ...
.. ..$tti: int [1:415, 1:3] 1 3 1 3 1 1 2 1 2 2 ... .. ..$ vv :Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
.. .. .. ..@ i       : int [1:1274] 26 37 61 95 192 207 20 30 55 127 ...
.. .. .. ..@ j       : int [1:1274] 0 0 0 0 0 0 1 1 1 1 ...
.. .. .. ..@ Dim     : int [1:2] 223 223
.. .. .. ..@ Dimnames:List of 2
.. .. .. .. ..$: NULL .. .. .. .. ..$ : NULL
.. .. .. ..@ x       : num [1:1274] 1 1 1 1 1 1 1 1 1 1 ...
.. .. .. ..@ factors : list()
..$segm :List of 2 .. ..$ bnd:List of 5
.. .. ..$loc : NULL .. .. ..$ idx   : int [1:29, 1:2] 193 218 211 217 194 195 196 197 221 213 ...
.. .. ..$grp : int [1:29, 1] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..$ is.bnd: logi TRUE
.. .. ..$crs : NULL .. .. ..- attr(*, "class")= chr "inla.mesh.segment" .. ..$ int:List of 5
.. .. ..$loc : NULL .. .. ..$ idx   : int [1:38, 1:2] 1 38 17 25 9 33 14 31 2 21 ...
.. .. ..$grp : int [1:38, 1] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..$ is.bnd: logi FALSE
.. .. ..$crs : NULL .. .. ..- attr(*, "class")= chr "inla.mesh.segment" ..$ idx     :List of 2
.. ..$loc : int [1:100] 39 40 41 42 43 44 45 46 47 48 ... .. ..$ lattice: NULL
..$crs : NULL ..- attr(*, "class")= chr "inla.mesh" - attr(*, "spde")=List of 7 ..$ model     : chr "matern"
..$n.spde : int 223 ..$ n.theta   : int 2
..$param.inla:List of 15 .. ..$ n            : int 223
.. ..$n.theta : int 2 .. ..$ M0           :Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
.. .. .. ..@ i       : int [1:223] 0 1 2 3 4 5 6 7 8 9 ...
.. .. .. ..@ j       : int [1:223] 0 1 2 3 4 5 6 7 8 9 ...
.. .. .. ..@ Dim     : int [1:2] 223 223
.. .. .. ..@ Dimnames:List of 2
.. .. .. .. ..$: NULL .. .. .. .. ..$ : NULL
.. .. .. ..@ x       : num [1:223] 190 171 208 117 129 ...
.. .. .. ..@ factors : list()
.. ..$M1 :Formal class 'dgTMatrix' [package "Matrix"] with 6 slots .. .. .. ..@ i : int [1:1497] 0 26 37 61 95 192 207 1 20 30 ... .. .. .. ..@ j : int [1:1497] 0 0 0 0 0 0 0 1 1 1 ... .. .. .. ..@ Dim : int [1:2] 223 223 .. .. .. ..@ Dimnames:List of 2 .. .. .. .. ..$ : NULL
.. .. .. .. ..$: NULL .. .. .. ..@ x : num [1:1497] 5.1604 -2.1409 -1.9881 -0.0821 -0.0883 ... .. .. .. ..@ factors : list() .. ..$ M2           :Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
.. .. .. ..@ i       : int [1:4041] 0 7 16 26 37 61 65 84 95 108 ...
.. .. .. ..@ j       : int [1:4041] 0 0 0 0 0 0 0 0 0 0 ...
.. .. .. ..@ Dim     : int [1:2] 223 223
.. .. .. ..@ Dimnames:List of 2
.. .. .. .. ..$: NULL .. .. .. .. ..$ : NULL
.. .. .. ..@ x       : num [1:4041] 0.1841 0.0137 0.0143 -0.1095 -0.1042 ...
.. .. .. ..@ factors : list()
.. ..$B0 : num [1:223, 1:3] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ B1           : num [1:223, 1:3] 0 0 0 0 0 0 0 0 0 0 ...
.. ..$B2 : num [1:223, 1:3] 1 1 1 1 1 1 1 1 1 1 ... .. .. ..- attr(*, "dimnames")=List of 2 .. .. .. ..$ : NULL
.. .. .. ..$: chr [1:3] "B" "" "" .. ..$ BLC          : num [1:6, 1:3] 0 0 0 0 -2.53 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$: chr [1:6] "theta.1" "theta.2" "tau.1" "kappa.1" ... .. .. .. ..$ : NULL
.. ..$theta.mu : num [1:2] 1.13 -2.4 .. ..$ theta.Q      : num [1:2, 1:2] 0.1 0 0 0.1
.. ..$transform : chr "identity" .. ..$ theta.initial: num [1:2] 1.13 -2.4
.. ..$fixed : logi [1:2] FALSE FALSE .. ..$ theta.fixed  : num(0)
..$f :List of 4 .. ..$ model          : chr "spde2"
.. ..$n : int 223 .. ..$ spde2.transform: chr "identity"
.. ..$hyper.default :List of 2 .. .. ..$ theta1:List of 3
.. .. .. ..$prior : chr "mvnorm" .. .. .. ..$ param  : num [1:6] 1.13 -2.4 0.1 0 0 ...
.. .. .. ..$initial: num 1.13 .. .. ..$ theta2:List of 1
.. .. .. ..$initial: num -2.4 ..$ BLC       : num [1:6, 1:3] 0 0 0 0 -2.53 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$: chr [1:6] "theta.1" "theta.2" "tau.1" "kappa.1" ... .. .. ..$ : NULL
..$mesh :List of 8 .. ..$ meta    :List of 5
.. .. ..$call : language inla.mesh.create(loc = rbind(loc, mesh2$loc), boundary = boundary[[2]],      interior = c(boundary2, interior2), extend = list(n = n[2], offset = offset[2]),  ...
.. .. ..$fmesher.args: chr "--input=input.s --cutoff=1e-12 --interior=input.segm.int.idx --interiorgrp=input.segm.int.grp --cet=16,-0.15 --rcdt=21,40,40 "| __truncated__ .. .. ..$ time        : num [1:5, 1:5] 0 0 0 0 0 0 0 0 0 0 ...
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$: chr [1:5] "pre" "fmesher" "post" "object" ... .. .. .. .. ..$ : NULL
.. .. ..$prefix : chr "/tmp/RtmpnLy1Nt/fmesher3e541a5fea57." .. .. ..$ is.refined  : logi TRUE
.. ..$manifold: chr "R2" .. ..$ n       : int 223
.. ..$loc : num [1:223, 1:3] 8.28 90.72 106.58 106.58 100.72 ... .. ..$ graph   :List of 5
.. .. ..$tv : int [1:415, 1:3] 29 127 73 64 200 26 50 49 12 23 ... .. .. ..$ vt : int [1:223, 1] 12 21 36 414 53 7 10 85 11 35 ...
.. .. ..$tt : int [1:415, 1:3] 26 212 161 91 53 372 69 113 65 111 ... .. .. ..$ tti: int [1:415, 1:3] 1 3 1 3 1 1 2 1 2 2 ...
.. .. ..$vv :Formal class 'dgTMatrix' [package "Matrix"] with 6 slots .. .. .. .. ..@ i : int [1:1274] 26 37 61 95 192 207 20 30 55 127 ... .. .. .. .. ..@ j : int [1:1274] 0 0 0 0 0 0 1 1 1 1 ... .. .. .. .. ..@ Dim : int [1:2] 223 223 .. .. .. .. ..@ Dimnames:List of 2 .. .. .. .. .. ..$ : NULL
.. .. .. .. .. ..$: NULL .. .. .. .. ..@ x : num [1:1274] 1 1 1 1 1 1 1 1 1 1 ... .. .. .. .. ..@ factors : list() .. ..$ segm    :List of 2
.. .. ..$bnd:List of 5 .. .. .. ..$ loc   : NULL
.. .. .. ..$idx : int [1:29, 1:2] 193 218 211 217 194 195 196 197 221 213 ... .. .. .. ..$ grp   : int [1:29, 1] 0 0 0 0 0 0 0 0 0 0 ...
.. .. .. ..$is.bnd: logi TRUE .. .. .. ..$ crs   : NULL
.. .. .. ..- attr(*, "class")= chr "inla.mesh.segment"
.. .. ..$int:List of 5 .. .. .. ..$ loc   : NULL
.. .. .. ..$idx : int [1:38, 1:2] 1 38 17 25 9 33 14 31 2 21 ... .. .. .. ..$ grp   : int [1:38, 1] 0 0 0 0 0 0 0 0 0 0 ...
.. .. .. ..$is.bnd: logi FALSE .. .. .. ..$ crs   : NULL
.. .. .. ..- attr(*, "class")= chr "inla.mesh.segment"
.. ..$idx :List of 2 .. .. ..$ loc    : int [1:100] 39 40 41 42 43 44 45 46 47 48 ...
.. .. ..$lattice: NULL .. ..$ crs     : NULL
.. ..- attr(*, "class")= chr "inla.mesh"
..- attr(*, "class")= chr [1:3] "inla.spde2" "inla.spde" "inla.model.class"
- attr(*, "Q")=Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. ..@ i       : int [1:4041] 0 7 16 26 37 61 65 84 95 108 ...
.. ..@ p       : int [1:224] 0 17 34 51 67 87 100 116 133 150 ...
.. ..@ Dim     : int [1:2] 223 223
.. ..@ Dimnames:List of 2
.. .. ..$: NULL .. .. ..$ : NULL
.. ..@ x       : num [1:4041] 7.43e+02 2.23e-05 2.33e-05 -3.41e-01 -3.17e-01 ...
.. ..@ factors : list()
- attr(*, "A")=Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. ..@ i       : int [1:300] 23 17 86 46 54 15 97 83 60 92 ...
.. ..@ p       : int [1:224] 0 1 2 2 3 3 3 3 5 7 ...
.. ..@ Dim     : int [1:2] 100 223
.. ..@ Dimnames:List of 2
.. .. ..$: NULL .. .. ..$ : NULL
.. ..@ x       : num [1:300] 0 0 0 0 0 0 0 0 0 0 ...
.. ..@ factors : list()
- attr(*, "cpu")=Class 'difftime'  atomic [1:3] 0.0691 0.0267 0.0958
.. ..- attr(*, "units")= Named chr "secs"
.. .. ..- attr(*, "names")= chr "prep"
## add squared exponential autoregressive structure
x <- x.k
for (j in 2:k) x[, j] <- rho * x[, j - 1] + sqrt(1 - rho^2) * x.k[, j]
dim(x)
[1] 100  10
library(dplyr)
library(reshape2)
x <- data.frame(x) %>% melt()
variable         value
1  sample1 -0.0006269303
2  sample1  0.0583550105
3  sample1  0.0205464030
4  sample1 -0.1838224593
5  sample1 -0.0550666120
6  sample1 -0.0263546934
## expand the coords into a data frame that is repeated for each time.
coords.st <- Reduce("rbind", lapply(1:k, function(x) data.frame(Year = x, coords)))
Year long lat
1    1   27  66
2    1   37  35
3    1   57  27
4    1   89  97
5    1   20  61
6    1   86  21
coords.st$y = x$value
Year long lat             y
1    1   27  66 -0.0006269303
2    1   37  35  0.0583550105
3    1   57  27  0.0205464030
4    1   89  97 -0.1838224593
5    1   20  61 -0.0550666120
6    1   86  21 -0.0263546934
ggplot(coords.st, aes(y = lat, x = long)) + geom_point(aes(color = y)) + facet_wrap(~Year)
spde <- attr(x.k, "spde")

### SPDE model

coords = as.data.frame(coords.st) %>% dplyr:::select(long, lat) %>% distinct %>%
as.matrix
mesh <- inla.mesh.2d(coords, max.edge = c(10, 40))
plot(mesh)
points(coords, col = "red", pch = 16)
spde <- inla.spde2.matern(mesh, alpha = 2)
s.index <- inla.spde.make.index("spatial.field", n.spde = spde$n.spde, n.group = k) tail(s.index$spatial.field.group)
[1] 10 10 10 10 10 10
tail(s.index$spatial.field.repl) [1] 1 1 1 1 1 1 str(s.index) List of 3$ spatial.field      : int [1:4520] 1 2 3 4 5 6 7 8 9 10 ...
$spatial.field.group: int [1:4520] 1 1 1 1 1 1 1 1 1 1 ...$ spatial.field.repl : int [1:4520] 1 1 1 1 1 1 1 1 1 1 ...
str(inla.spde.make.index("spatial.field", n.spde = spde$n.spde, n.repl = k)) List of 3$ spatial.field      : int [1:4520] 1 2 3 4 5 6 7 8 9 10 ...
$spatial.field.group: int [1:4520] 1 1 1 1 1 1 1 1 1 1 ...$ spatial.field.repl : int [1:4520] 1 1 1 1 1 1 1 1 1 1 ...
tail(inla.spde.make.index("spatial.field", n.spde = spde$n.spde, n.repl = k)$spatial.field.group)
[1] 1 1 1 1 1 1
tail(inla.spde.make.index("spatial.field", n.spde = spde$n.spde, n.repl = k)$spatial.field.repl)
[1] 10 10 10 10 10 10
## What is the difference between a group and a repl?

### The projector matrix

A <- inla.spde.make.A(mesh = mesh, loc = as.matrix(coords.st[, c("long", "lat")]),
group = coords.st$Year) ### The stack stack <- inla.stack(data = list(y = coords.st$y), A = list(A, 1), effects = list(s.index,
list(Intercept = rep(1, nrow(coords.st)))), tag = "est")

### Fit the INLA model

form <- y ~ 0 + Intercept + f(spatial.field, model = spde, group = spatial.field.group,
control.group = list(model = "ar1"))
data.inla <- inla(form, data = inla.stack.data(stack), control.predictor = list(compute = TRUE,
A = inla.stack.A(stack)), control.fixed = list(expand.factor.strategy = "inla"))
summary(data.inla)
Call:
c("inla(formula = form, data = inla.stack.data(stack), control.predictor = list(compute = TRUE, ",  "    A = inla.stack.A(stack)), control.fixed = list(expand.factor.strategy = \"inla\"))" )

Time used:
Pre-processing    Running inla Post-processing           Total
0.6240        565.6171          0.8941        567.1352

Fixed effects:
mean     sd 0.025quant 0.5quant 0.975quant   mode kld
Intercept 0.0034 0.0047    -0.0057   0.0034     0.0126 0.0034   0

Random effects:
Name	  Model
spatial.field   SPDE2 model

Model hyperparameters:
mean        sd 0.025quant
Precision for the Gaussian observations 39352.1888 2.395e+04 10581.3728
Theta1 for spatial.field                   -1.1877 2.387e+00    -4.7099
Theta2 for spatial.field                    1.0238 1.233e+00    -1.5756
GroupRho for spatial.field                  0.6844 2.810e-02     0.6305
0.5quant 0.975quant       mode
Precision for the Gaussian observations 33798.4024  1.013e+05 24399.5345
Theta1 for spatial.field                   -1.5211  3.773e+00    -3.9537
Theta2 for spatial.field                    1.1982  2.837e+00     2.4408
GroupRho for spatial.field                  0.6838  7.401e-01     0.6805

Expected number of effective parameters(std dev): 980.10(12.10)
Number of equivalent replicates : 1.02

Marginal log-Likelihood:  1531.66
Posterior marginals for linear predictor and fitted values computed
rf <- inla.spde2.result(data.inla, "spatial.field", spde, do.transf = TRUE)
# str(rf)

par(mfrow = c(2, 2), mar = c(3, 3, 1, 0.1), mgp = 2:0)
plot(data.inla$marginals.hyper[[3]], type = "l", xlab = "beta", ylab = "Density") abline(v = rho, col = 2) plot(rf$marginals.variance.nominal[[1]], type = "l", xlab = "sigma[x]", ylab = "Density")
abline(v = nu, col = 2)
plot(rf$marginals.kappa[[1]], type = "l", xlab = "kappa", ylab = "Density") abline(v = kappa, col = 2) plot(rf$marginals.range.nominal[[1]], type = "l", xlab = "range nominal", ylab = "Density")
abline(v = sqrt(8)/kappa, col = 2)

#### Explore the posterior of the random field

str(idat <- inla.stack.index(stack, "est")$data) int [1:1000] 1 2 3 4 5 6 7 8 9 10 ... cor(coords.st$y, data.inla$summary.linear.predictor$mean[idat])
[1] 0.9999621
stepsize <- 4 * 1/111
# nxy <- round(c(diff(range(as.vector(coords.st[,'long']))),
# diff(range(as.vector(coords.st[,'lat']))))/stepsize) projgrid <-
# inla.mesh.projector(mesh,
# xlim=range(as.vector(coords.st[,'long'])),ylim=range(as.vector(coords.st[,'lat'])),
# dims=nxy)
coords.sp <- SpatialPolygons(list(Polygons(list(Polygon(cbind(x = c(0, 100, 100,
0), y = c(0, 0, 100, 100)))), ID = "1")))
coords.sp <- spsample(coords.sp, n = 1000, type = "regular")
SpatialPoints:
x1        x2
[1,]  1.205305 0.1450983
[2,]  4.367583 0.1450983
[3,]  7.529861 0.1450983
[4,] 10.692138 0.1450983
[5,] 13.854416 0.1450983
[6,] 17.016694 0.1450983
Coordinate Reference System (CRS) arguments: NA
class(coords.sp)
[1] "SpatialPoints"
attr(,"package")
[1] "sp"
str(coords.sp)
Formal class 'SpatialPoints' [package "sp"] with 3 slots
..@ coords     : num [1:1024, 1:2] 1.21 4.37 7.53 10.69 13.85 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$: NULL .. .. ..$ : chr [1:2] "x1" "x2"
..@ bbox       : num [1:2, 1:2] 1.205 0.145 99.236 98.176
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$: chr [1:2] "x1" "x2" .. .. ..$ : chr [1:2] "min" "max"
..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
.. .. ..@ projargs: chr NA
projgrid <- inla.mesh.projector(mesh, loc = as.matrix(as.data.frame(coords.sp)))
## The prediction for each time can be done by
xmean <- list()
newdata <- NULL

for (j in 1:k) newdata = rbind(newdata, data.frame(Year = j, coords.sp, mean = inla.mesh.project(projgrid,
data.inla$summary.random$spatial.field$mean[s.index$spatial.field.group == j] +
data.inla$summary.fixed$mean[1])))

Year        x1        x2        mean
1    1  1.205305 0.1450983 0.003424414
2    1  4.367583 0.1450983 0.003424594
3    1  7.529861 0.1450983 0.003424886
4    1 10.692138 0.1450983 0.003425941
5    1 13.854416 0.1450983 0.003428833
6    1 17.016694 0.1450983 0.003431504
ggplot(newdata, aes(y = x2, x = x1)) + geom_tile(aes(fill = mean)) + facet_wrap(~Year)
Now lets fit including lat and long as random effects.
stack.e <- inla.stack(data = list(y = coords.st$y), A = list(A, 1, 1, 1), effects = list(s.index, list(Intercept = rep(1, nrow(coords.st))), list(long = coords.st$long), list(lat = coords.st$lat)), tag = "est") newdata = data.frame(long = rep(mesh$loc[, 1], k), lat = rep(mesh$loc[, 2], k), t = rep(1:k, each = mesh$n))

A.p = Diagonal(n = nrow(newdata))
stack.p = inla.stack(data = list(y = NA), A = list(A.p, 1, 1, 1), effects = list(s.index,
list(Intercept = rep(1, nrow(newdata))), list(lon = newdata$long), list(lat = newdata$lat)),
tag = "pred")
stack = inla.stack(stack.e, stack.p)
The fitted model then includes fixed effects of lat and long. Note also, that since we are going to want to access the estimated fitted values of the prediction set, we must include the compute=TRUE argument to the control.predictor() parameter.
formula = y ~ 0 + Intercept + lat + long + f(spatial.field, model = spde, group = spatial.field.group,
control.group = list(model = "ar1"))

data.inla <- inla(formula, data = inla.stack.data(stack), control.predictor = list(A = inla.stack.A(stack),
compute = TRUE))
summary(data.inla)
Call:
c("inla(formula = formula, data = inla.stack.data(stack), control.predictor = list(A = inla.stack.A(stack), ",  "    compute = TRUE))")

Time used:
Pre-processing    Running inla Post-processing           Total
0.7510        422.8323          1.5971        425.1803

Fixed effects:
mean     sd 0.025quant 0.5quant 0.975quant   mode kld
Intercept -8e-04 0.0113    -0.0230   -7e-04     0.0214 -7e-04   0
lat        1e-04 0.0001    -0.0002    1e-04     0.0004  1e-04   0
long       0e+00 0.0001    -0.0003    0e+00     0.0002  0e+00   0

Random effects:
Name	  Model
spatial.field   SPDE2 model

Model hyperparameters:
mean        sd 0.025quant
Precision for the Gaussian observations 39726.8987 2.456e+04 10712.2830
Theta1 for spatial.field                   -2.2313 8.455e-01    -4.0981
Theta2 for spatial.field                    1.5831 4.376e-01     0.8604
GroupRho for spatial.field                  0.6813 2.480e-02     0.6306
0.5quant 0.975quant       mode
Precision for the Gaussian observations 33893.0710  1.034e+05 24397.2518
Theta1 for spatial.field                   -2.1352 -8.379e-01    -1.7642
Theta2 for spatial.field                    1.5335  2.550e+00     1.3420
GroupRho for spatial.field                  0.6819  7.280e-01     0.6831

Expected number of effective parameters(std dev): 980.00(11.91)
Number of equivalent replicates : 1.02

Marginal log-Likelihood:  1507.90
Posterior marginals for linear predictor and fitted values computed
p.index = inla.stack.index(stack, "pred")$data data.mean = data.inla$summary.fitted.values$mean[p.index] data = as.data.frame(coords.st) coords.grid = expand.grid(long = seq(min(data$long), max(data$long), len = 100), lat = seq(min(data$lat), max(data$lat), len = 100)) data.inla.projector <- inla.mesh.projector(mesh, loc = as.matrix(coords.grid)) fit = sd = NULL for (i in 1:k) { ii = (i - 1) * mesh$n + 1
jj = i * mesh\$n
fit[[i]] = cbind(coords.grid, mean = c(inla.mesh.project(data.inla.projector,
data.mean[ii:jj])), Year = i)
}

fit = do.call("rbind", fit) %>% filter(!is.na(mean))