Fitting Occupancy models with R-INLA

This section describes the steps to fit occupancy models in R-INLA using simulated data (simulation details can be found in the Data Simulation tab).

Simple Spatial Occupancy Model

Model description goes in here …

Set up

We first load the data and prepare it in the format that is required by the INLA library.

Code

library(INLA)
library(inlabru)
library(fmesher)
library(tidyverse)
library(sf)
library(terra)
library(dplyr)

SSOM <- read.csv("Occ_data_1.csv")
x_covariate <- terra::rast('raster data/x_covariat.tif')
g_covariate <- terra::rast('raster data/g_covariat.tif')

# Extract the covariate values (NOTE: adapt this if inlabru is used)

# Convert to sf 
SSOM <- SSOM |>
  st_as_sf(coords = c('x.loc','y.loc')) 

#evaluate covariartes at each coordinate
SSOM = SSOM |> 
        mutate(terra::extract(x_covariate,st_coordinates(SSOM)),
               terra::extract(g_covariate,st_coordinates(SSOM)))

Table 1: First 6 entries of the occupancy data
cellid	y	nvisits	geometry	x_s	g_s
2	0	5	POINT (4.5 1.5)	0.0674760	1.4913634
5	3	4	POINT (13.5 1.5)	-0.2770668	0.9355086
6	2	3	POINT (16.5 1.5)	-0.4963150	0.8894699
7	4	5	POINT (19.5 1.5)	-0.4927090	0.7932032
9	0	1	POINT (25.5 1.5)	-0.2284233	0.5705246
23	0	4	POINT (67.5 1.5)	0.7247477	-0.1598709

Create the mesh … add details

Code

boundary_sf = st_bbox(c(xmin = 0, xmax = 300, ymax = 0, ymin = 300)) |> 
  st_as_sfc()

mesh = fm_mesh_2d(loc.domain = st_coordinates(boundary_sf)[,1:2],
                    offset = c(-0.1, -.2),
                    max.edge = c(15, 30))
matern <- inla.spde2.pcmatern(mesh,
                              prior.range = c(100, 0.5),  
                              prior.sigma = c(1, 0.5))

Create projector A matrix and make stacks .

add list of the arguments for building the stack (switch with inlabru details)

Code

# A_sp <- inla.spde.make.A(mesh = mesh, 
#                       loc = SSOM[,c('x.loc','y.loc')] |> as.matrix())

A_sp <- inla.spde.make.A(mesh = mesh, 
                      loc = st_coordinates(SSOM))

iset_sp <- inla.spde.make.index(name = "spatial_field", matern$n.spde)


stk <- inla.stack(data=list(Ycounts = SSOM$y, # observed occurrences
                            Ncounts = SSOM$nvisits, # number of visits
                            det_cov = SSOM$x_s, # detection covariate
                            Int_det = rep(1,length(SSOM$y))), # Baseline detection 
                  A=list(A_sp,1),  # the A matrix; the 1 is included to make the list(covariates)
                  effects=list(c(list(Int_occ=1), #the Intercept
                                 iset_sp),  #the spatial index
                               #the covariates
                               list(occ_cov = SSOM$x_s)), 
                  #this is a quick name so yo can call upon easily
                  tag='ssom')

Now we define the model components (left hand side -observational model components; right hand side - state process components) and fit the model (switch with inlabru details):

Code

formula_ssom <- inla.mdata(cbind(Ycounts,Ncounts),cbind(Int_det,det_cov)) ~  -1 + Int_occ +  occ_cov +  f(spatial_field, model=matern) 


model_ssom <- inla(formula_ssom, data=inla.stack.data(stk), 
                 family= '0binomialS',
                 control.fixed =  list(prec = 1/2.72, prec.intercept = 1/2.72),
                 control.predictor=list(A=inla.stack.A(stk),compute=TRUE),
                 control.compute = list(dic = TRUE, waic = TRUE, config = TRUE),
                 verbose = FALSE,
                 control.family = list(control.link = list(model = "logit"),
                                       link.simple = "logit",
                 hyper = list(
                   beta1 = list(param = c(0,1), initial = -1),
                   beta2 = list(param = c(0,1/2.72)))
                 )
                 )

Results

Show the summary results in Table 2:

Table 2: summary results from output
par	true	mean	quant0.025	quant0.975
\(\beta_0\)	-0.85	-2.34	-3.99	-1.00
\(\beta_1\)	1.50	1.14	0.81	1.47
\(\alpha_0\)	0.41	0.35	0.22	0.48
\(\alpha_1\)	1.00	0.11	-0.01	0.23
\(\rho\)	100.00	203.39	100.41	383.38
\(\sigma\)	1.00	1.55	0.98	2.40

show some plots: