Spatial Prediction with INLA

Introduction

This vignette illustrates the use of INLA for spatial prediction using examples from Blangiardo and Cameletti (2015) and a classic point process dataset. For prediction of continuous spatial processes, the Lindgren, Rue, and Lindström (2011) stochastic partial differential equations (SPDE) approach is used to approximate the process through an areal Gaussian Markov random field (GMRF) representation. Finally, Log-Gaussian Cox process models are fit using the pseudodata approach of Simpson et al. (2016).

GMRF Background

Blangiardo and Cameletti (2015) section 6.1.

• Observations aggregated to disjoint areal regions indexed by \(i\).
• Each region has unique parameter \(\theta_{i}\).
• \(\mathcal{N}(i)\) is the set of indices of neighbors of region \(i\) and \(\mathcal{N}_{i} = |\mathcal{N}(i)|\) is the number of neighbor of region \(i\).
• Local Markov property: given \(\boldsymbol{\theta}_{\mathcal{N}(i)}\), \(\theta_{i}\) is independent of all other \(\theta_{j}\).
• Then the precision matrix \(\mathbf{Q}\) of \(\boldsymbol{\theta}\) is sparse because only neighbors have nonzero coprecisions.
• Besag-York-Molliè model:
  • Exchangeable random effects \(u_{i}\) with intrinsic conditional autoregressive (iCAR) structure.
  • \(u_{i}|\mathbf{u}_{-i} \sim \mathrm{N}\left(\mu_{i} + \sum_{j} a_{ij} (u_{j} - \mu_{j}) / \mathcal{N}_{i}, \sigma_{u}^{2} / \mathcal{N}_{i}\right)\)
  • (What are the \(a_{ij}?\)).
  • iCAR is an improper prior because covariance matrix not positive definite but this is ok for random effects.

SPDE Background

• Spatial process \(Z(\mathbf{s})\) with \(d\)-dimensional domain is to be modeled.
• Solve \((\kappa^{2} - \Delta)^{\alpha / 2} (\tau Z(\mathbf{s})) = W(\mathbf{s})\).
  • \(\kappa\) is a range parameter.
  • \(\Delta\) is the Laplacian.
  • \(\alpha\) is a smoothness parameter.
  • \(\tau\) is a precision parameter.
  • \(W(\mathbf{s})\) is a Gaussian white noise process with mean 0 and variance 1.
  • Exact solution for stationary isotropic processes and integer \(\alpha\): \(Z(\mathbf{s})\) is a GMRF that well-approximates a Gaussian field with Matèrn covariance function with \(\nu = \alpha - d/2\).
    • Note for \(d = 2\) this excludes the exponential covariance (\(\nu = 1/2\)).
  • Some other results for non-stationary and anisotropic processes.
• Finite element approximation \(Z(\mathbf{s}) = \sum_{i = 1}^{n} z_{i} \psi_{i}(\mathbf{s})\) where \(\mathbf{z} = (z_{1}, \dots, z_{n})'\) is multivariate normal and \(\{\psi_{i}(\mathbf{s})\}_{i = 1}^{n}\) is a set of piecewise linear basis functions.
  • Specifically, \(\psi_{i}(\mathbf{s}) = 1\) at the \(i\)th node of a triangulation of the domain and 0 elsewhere.
  • \(Z(\mathbf{s})\) for arbitrary \(\mathbf{s}\) found by linear interpolation.
• Implemented as a black box in INLA.
  • User defines the triangulation but never sees the GMRF representation.
• Reduces covariance/precision matrix calculations from \(\mathcal{O}(n^{3})\) to \(\mathcal{O}(n^{3/2})\).

Geostatistics Example

Toy dataset from Blangiardo and Cameletti (2015). A spatial process is observed via a random sampling plan that concentrates observations in the lower left of the unit square.

# Plot the data.
plot(s2 ~ s1, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0), data = SPDEtoy, pch = 19, asp = 1, main = 'Toy Data')

# Create a mesh for the SPDE method and then plot it.
toy_mesh <- inla.mesh.2d(as.matrix(SPDEtoy[,c('s1', 's2')]), max.edge = c(0.1, 0.2))
plot(toy_mesh, asp = 1)
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)

# SPDE projector matrix for estimation.
A_est <- inla.spde.make.A(toy_mesh, as.matrix(SPDEtoy[,c('s1', 's2')]))

# Initialize exponential covariance structure for SPDE.
spde <- inla.spde2.matern(mesh = toy_mesh, alpha = 2)

# Set up stack for estimation.
stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = spde$n.spde)
stack_est <- inla.stack(data = list(y = SPDEtoy$y), A = list(A_est), effects = list(c(stack_index, list(intercept = 1))), tag = 'est')

# Create a grid for prediction.
toy_nx <- 50
toy_ny <- 50
toy_grid <- expand.grid(x = seq(0, 1, length.out = toy_nx), y = seq(0, 1, length.out = toy_ny))

# SPDE projector matrix for prediction.
A_pred <- inla.spde.make.A(mesh = toy_mesh, loc = as.matrix(toy_grid))

# Set up stacks for prediction.
stack_latent <- inla.stack(data = list(xi = NA), A = list(A_pred), effects = list(stack_index), tag = 'pred_latent')
stack_response <- inla.stack(data = list(y = NA), A = list(A_pred), effects = list(c(stack_index, list(intercept = 1))), tag = 'pred_response')

# Join all three stacks.
stacks <- inla.stack(stack_est, stack_latent, stack_response)

# Fit the model with INLA.
toy_fit <- inla(
  y ~ -1 + intercept + f(spatial_field, model = spde),
  data = inla.stack.data(stacks),
  control.predictor = list(A = inla.stack.A(stacks), compute = TRUE)
)

# Output posterior summaries.
toy_fit$summary.fixed

toy_fit$summary.hyperpar

# Extract posterior mean of latent spatial field.
index_latent <- inla.stack.index(stacks, tag = 'pred_latent')$data
post_mean <- toy_fit$summary.linear.predictor[index_latent, 'mean']
post_sd <- toy_fit$summary.linear.predictor[index_latent, 'sd']

# Plot the posterior mean and SD of the latent spatial field.
plot(im(matrix(post_mean, nrow = toy_nx, ncol = toy_ny), xrange = range(toy_grid$x), yrange = range(toy_grid$y)), main = 'Posterior Mean of Spatial Field')
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)

plot(im(matrix(post_sd, nrow = toy_nx, ncol = toy_ny), xrange = range(toy_grid$x), yrange = range(toy_grid$y)), main = 'Posterior SD of Spatial Field')
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)

Bei Dataset

Example from Møller and Waagepetersen (2007), Beilschmiedia pendula Lauraceae locations in a plot in Panama. bei dataset in spatstat (Baddeley and Turner 2005). The illustrations below fit a stationary LGCP model with no covariates, a Matèrn covariance function with \(\nu = 1\), and INLA’s default (flat?) priors for all parameters.

# Plot the full point pattern.
plot(bei, pch = '.', cols = 'black', main = 'Realized Point Pattern')

# Take a sample of quadrats and plot the observed point pattern.
set.seed(84323)
N_QUADS <- 10
QUAD_SIZE <- 50

w_edge <- Frame(bei)$xrange[1]
e_edge <- Frame(bei)$xrange[2]
s_edge <- Frame(bei)$yrange[1]
n_edge <- Frame(bei)$yrange[2]

botleft <- cbind(
  runif(N_QUADS, w_edge, e_edge - QUAD_SIZE),
  runif(N_QUADS, s_edge, n_edge - QUAD_SIZE)
)
bei_interior <- lapply(seq_len(nrow(botleft)), function(r){return(
    cbind(
      botleft[r, 1] + c(0, 0, QUAD_SIZE, QUAD_SIZE),
      botleft[r, 2] + c(0, QUAD_SIZE, QUAD_SIZE, 0)
    )
  )})
bei_win <- do.call(
  union.owin,
  apply(botleft, 1, function(x){return(
    owin(x[1] + c(0, QUAD_SIZE), x[2] + c(0, QUAD_SIZE))
  )})
)
bei_hole <- bei[complement.owin(bei_win, frame = Frame(bei))]
bei_samp <- bei[bei_win]
bei_window_full <- Window(bei)

plot(bei_hole, main = 'Observed Region with Holes', pch = '.', cols = 'black')

plot(bei_window_full, main = 'Observed Sample')
plot(bei_win, add = TRUE)
plot(bei_samp, pch = '.', cols = 'black', add = TRUE)

Triangulation Meshes for SPDE Approach

The mesh should be fine in observed regions to accurately represent complex local structure but can be coarsened in unobserved regions where the model cannot infer as much detail. I am including margins to explore model behavior away from the observed regions.

# Parameters to experiment with.
MAX_EDGE_LENGTH <- 25
MAX_EDGE_EXT <- 50
MARGIN <- 100

# Mesh covering the site.
bei_boundary <- inla.mesh.segment(loc = do.call(cbind, vertices.owin(Window(bei))))
bei_full_mesh <- inla.mesh.create(
  boundary = bei_boundary,
  refine = list(max.edge = MAX_EDGE_LENGTH)
)
bei_full_spde <- inla.spde2.matern(mesh = bei_full_mesh)
plot(bei_full_mesh, asp = 1, main = 'Fine Mesh')
points(bei, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh')

# Mesh including a margin outside the site.
margin_mesh <- inla.mesh.2d(
  loc = bei_full_mesh$loc[,1:2], # Include nodes from site.
  offset = MARGIN,
  max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
margin_spde <- inla.spde2.matern(mesh = margin_mesh)
plot(margin_mesh, asp = 1, main = 'Fine Mesh with Coarse Margin')
points(bei, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh with Coarse Margin')

# Meshs with coarser resolution in quadrats.
quad_hole <- do.call(
  inla.mesh.segment,
  lapply(seq_along(bei_interior), function(i){
    return(inla.mesh.segment(loc = bei_interior[[i]], grp = i - 1))
  })
)
bei_hole_mesh0 <- inla.mesh.create(
  boundary = list(bei_boundary, quad_hole),
  refine = list(max.edge = MAX_EDGE_LENGTH)
)
bei_hole0_spde <- inla.spde2.matern(mesh = bei_hole_mesh0)
plot(bei_hole_mesh0, asp = 1, main = 'Fine Mesh with Holes')
points(bei_hole, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh with Holes')

bei_hole_mesh <- inla.mesh.create(
  loc = bei_hole_mesh0$loc[,1:2], # Include nodes from mesh with holes.
  boundary = bei_boundary,
  refine = list(max.edge = MAX_EDGE_EXT) # Fill in the rest with a coarser triangulation.
)
bei_hole_spde <- inla.spde2.matern(mesh = bei_hole_mesh)
plot(bei_hole_mesh, asp = 1, main = 'Fine Mesh with Coarse Holes')
points(bei_hole, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh with Coarse Holes')

# Meshs with finer resolution in quadrats.
quad_bnd <- do.call(
  inla.mesh.segment,
  lapply(seq_along(bei_interior), function(i){
    return(inla.mesh.segment(loc = apply(bei_interior[[i]], 2, rev), grp = i - 1))
  })
)
bei_samp_mesh0 <- inla.mesh.create(
  boundary = quad_bnd,
  refine = list(max.edge = MAX_EDGE_LENGTH)
)
bei_samp0_spde <- inla.spde2.matern(mesh = bei_samp_mesh0)
plot(bei_samp_mesh0, asp = 1, main = 'Fine Mesh in Quadrats')
points(bei_samp, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh in Quadrats')

bei_samp_mesh <- inla.mesh.create(
  loc = bei_samp_mesh0$loc[,1:2], # Include nodes from mesh in quads.
  boundary = bei_boundary,
  refine = list(max.edge = MAX_EDGE_EXT) # Fill in the rest with a coarser triangulation.
)
bei_samp_spde <- inla.spde2.matern(mesh = bei_samp_mesh)
plot(bei_samp_mesh, asp = 1, main = 'Coarse Mesh with Fine Quadrats')
points(bei_samp, pch = 19, cex = 0.25, col = 'red')
title('Coarse Mesh with Fine Quadrats')

# Meshes with varying resolution in quadrats and a margin.
margin_hole <- inla.mesh.2d(
  loc = bei_hole_mesh$loc[,1:2], # Include nodes from mesh with holes.
  offset = MARGIN,
  max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
margin_hole_spde <- inla.spde2.matern(mesh = margin_hole)
plot(margin_hole, asp = 1, main = 'Coarse Holes and Margin')
points(bei_hole, pch = 19, cex = 0.25, col = 'red')
title('Coarse Holes and Margin')

margin_samp <- inla.mesh.2d(
  loc = bei_samp_mesh$loc[,1:2], # Include nodes from quads.
  offset = MARGIN,
  max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
margin_samp_spde <- inla.spde2.matern(mesh = margin_samp)
plot(margin_samp, asp = 1, main = 'Fine Quadrats and Coarse Margin')
points(bei_samp, pch = 19, cex = 0.25, col = 'red')
title('Fine Quadrats and Coarse Margin')

# Identify which mesh nodes are in the oberserved region
# and create SPDE projectors for spatial mapping.
NPIX_X <- 400
NPIX_Y <- 200

obs_full <- rep(0, margin_mesh$n)
obs_full[inla.over_sp_mesh(as(Window(bei), 'SpatialPolygons'), margin_mesh, 'vertex')] <- 1
proj_margin_mesh <- inla.mesh.projector(margin_mesh, dims = c(NPIX_X, NPIX_Y))

obs_hole <- rep(0, margin_hole$n)
obs_hole[inla.over_sp_mesh(as(Window(bei_hole), 'SpatialPolygons'), margin_hole, 'vertex')] <- 1
proj_margin_hole <- inla.mesh.projector(margin_hole, dims = c(NPIX_X, NPIX_Y))

obs_samp <- rep(0, margin_samp$n)
obs_samp[inla.over_sp_mesh(as(Window(bei_samp), 'SpatialPolygons'), margin_samp, 'vertex')] <- 1
proj_margin_samp <- inla.mesh.projector(margin_samp, dims = c(NPIX_X, NPIX_Y))

Bei Dataset with gridding

Count the points in grid cells and fit a Poisson GLMM with the observed area in the cell as an exposure variable, as done in Illian, Sørbye, and Rue (2012) and many other examples. The plots of the cell counts are blank (white) where cells have no observed area.

NGRID_X <- 40
NGRID_Y <- 20

centers <- gridcenters(
  dilation(bei_window_full, max(NGRID_X, NGRID_Y)),
  NGRID_X, NGRID_Y)
dx <- sum(unique(centers$x)[1:2] * c(-1, 1)) / 2
dy <- sum(unique(centers$y)[1:2] * c(-1, 1)) / 2

bei_df <- data.frame(x = centers$x, y = centers$y,
                     count = NA_integer_, area = NA_real_)

system.time(
for(r in seq_len(nrow(bei_df))){
  bei_df$count[r] <- sum(bei$x >= bei_df$x[r] - dx &
                         bei$x < bei_df$x[r] + dx &
                         bei$y >= bei_df$y[r] - dy &
                         bei$y < bei_df$y[r] + dy)
  bei_df$area[r] <- area(Window(bei)[owin(c(bei_df$x[r] - dx, bei_df$x[r] + dx), c(bei_df$y[r] - dy, bei_df$y[r] + dy))])
}
)

   user  system elapsed 
  0.472   0.004   0.474 

par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(ifelse(bei_df$area > 0, bei_df$count, NA), nrow = length(unique(bei_df$x)))), unique(bei_df$x), unique(bei_df$y), unitname = 'meters'), ncolcours = range(bei_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(bei_window_full, border = 'white', add = TRUE)
points(bei, pch = '.', col = 'black')

# SPDE projector matrix for estimation.
full_A_est <- inla.spde.make.A(
  margin_mesh,
  as.matrix(bei_df[bei_df$area > 0, c('x', 'y')])
)

# Set up stack for estimation.
full_stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = margin_spde$n.spde)
full_stack_est <- inla.stack(data = list(count = bei_df$count[bei_df$area > 0], larea = log(bei_df$area[bei_df$area > 0])), A = list(full_A_est), effects = list(c(full_stack_index, list(intercept = 1))), tag = 'est')

# SPDE projector matrix for prediction.
full_A_pred <- inla.spde.make.A(mesh = margin_mesh, loc = as.matrix(bei_df[,c('x', 'y')]))

# Set up stacks for prediction.
full_stack_latent <- inla.stack(data = list(xi = NA), A = list(full_A_pred), effects = list(full_stack_index), tag = 'pred_latent')
full_stack_response <- inla.stack(data = list(count = NA), A = list(full_A_pred), effects = list(c(full_stack_index, list(intercept = 1))), tag = 'pred_response')

# Join all three stacks.
full_stacks <- inla.stack(full_stack_est, full_stack_latent, full_stack_response)

# Fit the model with INLA.
system.time(
bei_full_fit <- inla(
  count ~ -1 + intercept + f(spatial_field, model = margin_spde),
  offset = larea, family = 'poisson',
  data = inla.stack.data(full_stacks),
  control.predictor = list(A = inla.stack.A(full_stacks), compute = TRUE),
  verbose = TRUE
)
)

   user  system elapsed 
904.840   0.452 122.519 

# Output posterior summaries.
bei_full_fit$summary.fixed

bei_full_fit$summary.hyperpar

# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_mesh, bei_full_fit$summary.fixed$mean + bei_full_fit$summary.random$spatial_field$mean)),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior Mean of Log-Intensity')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_mesh, sqrt(bei_full_fit$summary.fixed$sd^2 + bei_full_fit$summary.random$spatial_field$sd^2))),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior SD of Log-Intensity')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

bei_hole_df <- data.frame(x = centers$x, y = centers$y,
                          count = NA_integer_, area = NA_real_)

system.time(
for(r in seq_len(nrow(bei_hole_df))){
  bei_hole_df$count[r] <- sum(bei_hole$x >= bei_hole_df$x[r] - dx &
                              bei_hole$x < bei_hole_df$x[r] + dx &
                              bei_hole$y >= bei_hole_df$y[r] - dy &
                              bei_hole$y < bei_hole_df$y[r] + dy)
  bei_hole_df$area[r] <- area(Window(bei_hole)[owin(c(bei_hole_df$x[r] - dx, bei_hole_df$x[r] + dx), c(bei_hole_df$y[r] - dy, bei_hole_df$y[r] + dy))])
}
)

   user  system elapsed 
  1.168   0.000   1.170 

par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(ifelse(bei_hole_df$area > 0, bei_hole_df$count, NA), nrow = length(unique(bei_hole_df$x)))), unique(bei_hole_df$x), unique(bei_hole_df$y), unitname = 'meters'), ncolcours = range(bei_hole_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = '#00000040')

# SPDE projector matrix for estimation.
hole_A_est <- inla.spde.make.A(
  margin_hole,
  as.matrix(bei_hole_df[bei_hole_df$area > 0, c('x', 'y')])
)

# Set up stack for estimation.
hole_stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = margin_hole_spde$n.spde)
hole_stack_est <- inla.stack(data = list(count = bei_hole_df$count[bei_hole_df$area > 0], larea = log(bei_hole_df$area[bei_hole_df$area > 0])), A = list(hole_A_est), effects = list(c(hole_stack_index, list(intercept = 1))), tag = 'est')

# SPDE projector matrix for prediction.
hole_A_pred <- inla.spde.make.A(mesh = margin_hole, loc = as.matrix(bei_hole_df[,c('x', 'y')]))

# Set up stacks for prediction.
hole_stack_latent <- inla.stack(data = list(xi = NA), A = list(hole_A_pred), effects = list(hole_stack_index), tag = 'pred_latent')
hole_stack_response <- inla.stack(data = list(count = NA), A = list(hole_A_pred), effects = list(c(hole_stack_index, list(intercept = 1))), tag = 'pred_response')

# Join all three stacks.
hole_stacks <- inla.stack(hole_stack_est, hole_stack_latent, hole_stack_response)

# Fit the model with INLA.
system.time(
bei_hole_fit <- inla(
  count ~ -1 + intercept + f(spatial_field, model = margin_hole_spde),
  offset = larea, family = 'poisson',
  data = inla.stack.data(hole_stacks),
  control.predictor = list(A = inla.stack.A(hole_stacks), compute = TRUE),
  verbose = TRUE
)
)

   user  system elapsed 
833.072   0.320 113.165 

# Output posterior summaries.
bei_hole_fit$summary.fixed

bei_hole_fit$summary.hyperpar

# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_hole, bei_hole_fit$summary.fixed$mean + bei_hole_fit$summary.random$spatial_field$mean)),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior Mean of Log-Intensity')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_hole, sqrt(bei_hole_fit$summary.fixed$sd^2 + bei_hole_fit$summary.random$spatial_field$sd^2))),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior SD of Log-Intensity')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

bei_samp_df <- data.frame(x = centers$x, y = centers$y,
                          count = NA_integer_, area = NA_real_)

system.time(
for(r in seq_len(nrow(bei_samp_df))){
  bei_samp_df$count[r] <- sum(bei_samp$x >= bei_samp_df$x[r] - dx &
                              bei_samp$x < bei_samp_df$x[r] + dx &
                              bei_samp$y >= bei_samp_df$y[r] - dy &
                              bei_samp$y < bei_samp_df$y[r] + dy)
  bei_samp_df$area[r] <- area(Window(bei_samp)[owin(c(bei_samp_df$x[r] - dx, bei_samp_df$x[r] + dx), c(bei_samp_df$y[r] - dy, bei_samp_df$y[r] + dy))])
}
)

   user  system elapsed 
  1.060   0.000   1.059 

par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(ifelse(bei_samp_df$area > 0, bei_samp_df$count, NA), nrow = length(unique(bei_samp_df$x)))), unique(bei_samp_df$x), unique(bei_samp_df$y), unitname = 'meters'), ncolcours = range(bei_samp_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = '#00000040')

# SPDE projector matrix for estimation.
samp_A_est <- inla.spde.make.A(
  margin_samp,
  as.matrix(bei_samp_df[bei_samp_df$area > 0, c('x', 'y')])
)

# Set up stack for estimation.
samp_stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = margin_samp_spde$n.spde)
samp_stack_est <- inla.stack(data = list(count = bei_samp_df$count[bei_samp_df$area > 0], larea = log(bei_samp_df$area[bei_samp_df$area > 0])), A = list(samp_A_est), effects = list(c(samp_stack_index, list(intercept = 1))), tag = 'est')

# SPDE projector matrix for prediction.
samp_A_pred <- inla.spde.make.A(mesh = margin_samp, loc = as.matrix(bei_samp_df[,c('x', 'y')]))

# Set up stacks for prediction.
samp_stack_latent <- inla.stack(data = list(xi = NA), A = list(samp_A_pred), effects = list(samp_stack_index), tag = 'pred_latent')
samp_stack_response <- inla.stack(data = list(count = NA), A = list(samp_A_pred), effects = list(c(samp_stack_index, list(intercept = 1))), tag = 'pred_response')

# Join all three stacks.
samp_stacks <- inla.stack(samp_stack_est, samp_stack_latent, samp_stack_response)

# Fit the model with INLA.
system.time(
bei_samp_fit <- inla(
  count ~ -1 + intercept + f(spatial_field, model = margin_samp_spde),
  offset = larea, family = 'poisson',
  data = inla.stack.data(samp_stacks),
  control.predictor = list(A = inla.stack.A(samp_stacks), compute = TRUE),
  verbose = TRUE
)
)

   user  system elapsed 
275.936   0.284  38.750 

# Output posterior summaries.
bei_samp_fit$summary.fixed

bei_samp_fit$summary.hyperpar

# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_samp, bei_samp_fit$summary.fixed$mean + bei_samp_fit$summary.random$spatial_field$mean)),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior Mean of Log-Intensity')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_samp, sqrt(bei_samp_fit$summary.fixed$sd^2 + bei_samp_fit$summary.random$spatial_field$sd^2))),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior SD of Log-Intensity')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

Bei Dataset with Simpson et al. (2016) method

This method relies upon the Lindgren, Rue, and Lindström (2011) approximation of the latent Gaussian field as a linear combination of a finite number of basis functions represented as a GMRF on the nodes of a triangulation of the space. Simpson et al. (2016) use the triangulation for numerical integration of the intensity function and show that the LGCP likelihood factors into the joint distribution of independent Poisson random variables corresponding to the points of the point pattern and the nodes of the triangulation. The model fitting proceeds using INLA to fit a Poisson model to pseudodata.

The pseudodata are constructed as follows.

• Let \(n\) be the size of the point pattern.
• Let \(p\) be the number of nodes of triangulation.
• Define the pseudodata as \(\mathbf{y} = (y_{1}, \dots, y_{p}, y_{p+1}, \dots, y_{p+n})'\) where \(y_{i} = 0\) for \(i = 1, \dots, p\) (the traingulation nodes) and \(y_{i} = 1\) for \(i = p+1, \dots, p+n\) (the observed points).
• Define \(\boldsymbol{\alpha} = (\alpha_{i}, \dots, \alpha_{p}, \alpha_{p+1}, \dots, \alpha_{p+n})'\) to encode the numerical integration scheme where \(\alpha_{i}\) is the numerical integration weight for \(i = 1, \dots, p\) (the traingulation nodes) and \(\alpha_{i} = 0\) for \(i = p+1, \dots, p+n\) (the observed points).

Then \(y_{i} \sim Poisson(\alpha_{i}\eta_{i})\) where \(\log(\eta_{i})\) is the SPDE representation of the GF at the location of the \(i\)th pseudodatum. See the paper for tedious notation regarding the definition of \(\eta_{i}\). Ultimately, the nodes become Poisson random variables with means equal to their respective integration weights times the intensity at that their locations (so the integration weight is an exposure variable), observed points become Poisson random variables with means of 1, and the likelihood is approximately proportional to

\[\prod_{i=1}^{p+n} \eta_{i}^{y_{i}} \exp(-\alpha_{i} \eta_{i}).\]

(Is there a missing \(\alpha_{i}\)?)

Covariates can easily be included when they are observed at the mesh nodes. Known “sampling effort” is accomodated by scaling the integration weights by the probability that a point would have been observed at that node (given the sampling plan), i.e. nodes in observed regions have the \(\alpha_{i}\) defined above and nodes in unobserved regions have \(\alpha_{i} = 0\). Weights can be scaled by values other than 0 or 1 to account for e.g. distance sampling with a (known) detection function or a (known) false negative rate for detection equipment. This adjustment assumes the point process of interest is observed through a known thinning process and then allows inference back to the intensity function of the unthinned process. More complicated detection processes are possible, e.g. Yuan et al. (2017) fit an LGCP to data obtianed through distance sampling while using splines to model the detection function.

full_pts <- cbind(bei$x, bei$y)

# Contruct the SPDE A matrix for nodes and points.
full_nV <- margin_mesh$n
full_nData <- dim(full_pts)[1]
full_LocationMatrix <- inla.mesh.project(margin_mesh, full_pts)$A
full_IntegrationMatrix <- sparseMatrix(i = 1:full_nV, j = 1:full_nV, x = rep(1, full_nV))
full_ObservationMatrix <- rbind(full_IntegrationMatrix, full_LocationMatrix)

# Get the integration weights.
full_IntegrationWeights <- diag(inla.mesh.fem(margin_mesh)$c0)
full_E_point_process <- c(obs_full * full_IntegrationWeights, rep(0, full_nData))

# Create the psuedodata.
full_fake_data <- c(rep(0, full_nV), rep(1, full_nData))

# Fit model to full site.
full_formula <- y ~ -1 + intercept + f(idx, model = margin_spde) # No covariates.
full_data <- list(y = full_fake_data, idx = 1:full_nV, intercept = rep(1, full_nV))

system.time(
result_full <- inla(
  formula = full_formula,
  data = full_data,
  family = 'poisson',
  control.predictor = list(A = full_ObservationMatrix),
  E = full_E_point_process,
  verbose = TRUE
)
)

   user  system elapsed 
446.200   0.532  65.662 

result_full$summary.fixed

result_full$summary.hyperpar

# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_mesh, result_full$summary.fixed$mean + result_full$summary.random$idx$mean)),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior Mean of Log-Intensity')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_mesh, sqrt(result_full$summary.fixed$sd^2 + result_full$summary.random$idx$sd^2))),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior SD of Log-Intensity')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

hole_pts <- cbind(bei_hole$x, bei_hole$y)

# Contruct the SPDE A matrix for nodes and points.
hole_nV <- margin_hole$n
hole_nData <- dim(hole_pts)[1]
hole_LocationMatrix <- inla.mesh.project(margin_hole, hole_pts)$A
hole_IntegrationMatrix <- sparseMatrix(i = 1:hole_nV, j = 1:hole_nV, x = rep(1, hole_nV))
hole_ObservationMatrix <- rbind(hole_IntegrationMatrix, hole_LocationMatrix)

# Get the integration weights.
hole_IntegrationWeights <- diag(inla.mesh.fem(margin_hole)$c0)
hole_E_point_process <- c(obs_hole * hole_IntegrationWeights, rep(0, hole_nData))

# Create the psuedodata.
hole_fake_data <- c(rep(0, hole_nV), rep(1, hole_nData))

# Fit model to site with holes.
hole_formula <- y ~ -1 + intercept + f(idx, model = margin_hole_spde) # No covariates.
hole_data <- list(y = hole_fake_data, idx = 1:hole_nV, intercept = rep(1, hole_nV))

system.time(
result_hole <- inla(
  formula = hole_formula,
  data = hole_data,
  family = 'poisson',
  control.predictor = list(A = hole_ObservationMatrix),
  E = hole_E_point_process,
  verbose = TRUE
)
)

   user  system elapsed 
383.692   0.448  56.575 

result_hole$summary.fixed

result_hole$summary.hyperpar

# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_hole, result_hole$summary.fixed$mean + result_hole$summary.random$idx$mean)),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior Mean of Log-Intensity')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_hole, sqrt(result_hole$summary.fixed$sd^2 + result_hole$summary.random$idx$sd^2))),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior SD of Log-Intensity')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

samp_pts <- cbind(bei_samp$x, bei_samp$y)

# Contruct the SPDE A matrix for nodes and points.
samp_nV <- margin_samp$n
samp_nData <- dim(samp_pts)[1]
samp_LocationMatrix <- inla.mesh.project(margin_samp, samp_pts)$A
samp_IntegrationMatrix <- sparseMatrix(i = 1:samp_nV, j = 1:samp_nV, x = rep(1, samp_nV))
samp_ObservationMatrix <- rbind(samp_IntegrationMatrix, samp_LocationMatrix)

# Get the integration weights.
samp_IntegrationWeights <- diag(inla.mesh.fem(margin_samp)$c0)
samp_E_point_process <- c(obs_samp * samp_IntegrationWeights, rep(0, samp_nData))

# Create the psuedodata.
samp_fake_data <- c(rep(0, samp_nV), rep(1, samp_nData))

# Fit model to quadrat-sampled site.
samp_formula <- y ~ -1 + intercept + f(idx, model = margin_samp_spde) # No covariates.
samp_data <- list(y = samp_fake_data, idx = 1:samp_nV, intercept = rep(1, samp_nV))

system.time(
result_samp <- inla(
  formula = samp_formula,
  data = samp_data,
  family = 'poisson',
  control.predictor = list(A = samp_ObservationMatrix),
  E = samp_E_point_process,
  verbose = TRUE
)
)

   user  system elapsed 
348.300   0.808  52.568 

result_samp$summary.fixed

result_samp$summary.hyperpar

# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_samp, result_samp$summary.fixed$mean + result_samp$summary.random$idx$mean)),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior Mean of Log-Intensity')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_samp, sqrt(result_samp$summary.fixed$sd^2 + result_samp$summary.random$idx$sd^2))),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior SD of Log-Intensity')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

Bei Dataset and inlabru

inlabru has a wrapper function to fit an LGCP with INLA, with a relatively easy-to-use interface for defining models and predicting arbitrary functions of latent variables. However, it is poorly documented, slow, and the documentation seems to imply that it does not support variable sampling effort (even though this appears to work).

bei_full_spdf <- as.SpatialPoints.ppp(bei)
cmp_full <- coordinates ~ mySmooth(map = coordinates, model = margin_spde) + Intercept

system.time(
bei_full_lgcp <- lgcp(cmp_full, bei_full_spdf, E = obs_full, options = list(verbose = TRUE))
)

    user   system  elapsed 
3776.056    1.712  504.862 

bei_full_lgcp$summary.fixed

bei_full_lgcp$summary.hyperpar

# Plot posterior means and posterior sd.
lambda_full <- predict(bei_full_lgcp, pixels(margin_mesh), ~ Intercept + mySmooth)
plot(lambda_full, main = 'Posterior Mean of Log-Intensity')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

plot(lambda_full['sd'] ,main = 'Posterior SD of Log-Intensity')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

bei_hole_spdf <- as.SpatialPoints.ppp(bei_hole)
cmp_hole <- coordinates ~ mySmooth(map = coordinates, model = margin_hole_spde) + Intercept

system.time(
bei_hole_lgcp <- lgcp(cmp_hole, bei_hole_spdf, E = obs_hole, options = list(verbose = TRUE))
)

    user   system  elapsed 
3369.264    1.652  437.205 

bei_hole_lgcp$summary.fixed

bei_hole_lgcp$summary.hyperpar

# Plot posterior means and posterior sd.
lambda_hole <- predict(bei_hole_lgcp, pixels(margin_hole), ~ Intercept + mySmooth)
plot(lambda_hole, main = 'Posterior Mean of Log-Intensity')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

plot(lambda_hole['sd'], main = 'Posterior SD of Log-Intensity')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

bei_samp_spdf <- as.SpatialPoints.ppp(bei_samp)
cmp_samp <- coordinates ~ mySmooth(map = coordinates, model = margin_samp_spde) + Intercept

system.time(
bei_samp_lgcp <- lgcp(cmp_samp, bei_samp_spdf, E = obs_samp, options = list(verbose = TRUE))
)

   user  system elapsed 
689.148   1.480 101.756 

bei_samp_lgcp$summary.fixed

bei_samp_lgcp$summary.hyperpar

# Plot posterior means and posterior sd.
lambda_samp <- predict(bei_samp_lgcp, pixels(margin_samp), ~ Intercept + mySmooth)
plot(lambda_samp, main = 'Posterior Mean of Log-Intensity')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

plot(lambda_samp['sd'], main = 'Posterior SD of Log-Intensity')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

Summary

All three methods give similar results for the full dataset and the dataset with holes, even when the gridding method uses a coarse grid. The intercept and random effects are shifted from method to method, but these are not separately indentifiable and the shifts cancel each other out. The methods each have different artifacts and edge effects apparent in the results from the sampled dataset. The pseudodata approach is the fastest except when a very coarse grid is used and a small region is observed.

References

Baddeley, Adrian, and Rolf Turner. 2005. “Spatstat: An R Package for Analyzing Spatial Point Patterns.” Journal of Statistical Software 12 (6): 1–42.

Blangiardo, Marta, and Michela Cameletti. 2015. Spatial and Spatio-Temporal Bayesian Models with R-INLA. Wiley.

Illian, Janine B, Sigrunn H Sørbye, and Håvard Rue. 2012. “A Toolbox for Fitting Complex Spatial Point Process Models Using Integrated Nested Laplace Approximation (Inla).” The Annals of Applied Statistics, 1499–1530.

Lindgren, Finn, Håvard Rue, and Johan Lindström. 2011. “An Explicit Link Between Gaussian Fields and Gaussian Markov Random Fields: The Stochastic Partial Differential Equation Approach.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73 (4): 423–98.

Møller, J, and RP Waagepetersen. 2007. “Modern Spatial Point Process Modelling and Inference.” Scandinavian Journal of Statistics 34: 643–711.

Simpson, Daniel, Janine B Illian, Finn Lindgren, Sigrunn H Sørbye, and Havard Rue. 2016. “Going Off Grid: Computationally Efficient Inference for Log-Gaussian Cox Processes.” Biometrika 103 (1): 49–70.

Yuan, Yuan, Fabian E Bachl, Finn Lindgren, David L Borchers, Janine B Illian, Stephen T Buckland, Haavard Rue, Tim Gerrodette, and others. 2017. “Point Process Models for Spatio-Temporal Distance Sampling Data from a Large-Scale Survey of Blue Whales.” The Annals of Applied Statistics 11 (4): 2270–97.