Spatiotemporal example

Studying the complex ecological systems across both space and time is imperative in understanding the full dynamics of species. This vignette illustrates the construction of a spatiotemporal ISDM using PointedSDMs, using data of species Colinus virginianus across Alabama (United States of America). The first step in this vignette is to load the required packages:


library(PointedSDMs)
library(inlabru)
library(ggplot2)
library(spocc)
library(INLA)
library(dplyr)
library(sp)
library(sf)


proj <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

AL <- USAboundaries::us_states(states = "Alabama", resolution = 'high')
AL <- as(AL, "sf")
st_crs(AL) <- proj

mesh <- inla.mesh.2d(boundary = inla.sp2segment(AL[1]), 
                     cutoff = 0.1,
                     max.edge = c(0.2, 0.8), 
                     offset = c(0.1, 0.2),
                     crs = st_crs(proj))

The first dataset we consider is obtained from the North American Breeding Bird Survey. This dataset may be loaded directly from the package, and contains observations of the species between 2015 and 2017. This dataset is treated as replicate present-absent, where every point is assumed to be a visited site (or route).


data("BBSColinusVirginianus")

The second dataset considered is obtained via the citizen science program, eBird. These data are obtained via the R package, spocc using the script below, where a separate object of data points was created for each year to ensure that the number of records per year is equal.


eBird2015 <- spocc::occ(
  query = 'Colinus virginianus',
  from = 'gbif',
  date = c("2015-01-01", "2015-12-31"),
  geometry = AL
)$gbif

eBird2016 <- spocc::occ(
  query = 'Colinus virginianus',
  from = 'gbif',
  date = c("2016-01-01", "2016-12-31"),
  geometry = AL
)$gbif

eBird2017 <- spocc::occ(
  query = 'Colinus virginianus',
  from = 'gbif',
  date = c("2017-01-01", "2017-12-31"),
  geometry = AL
)$gbif

eBird <- data.frame(eBird2015$data[[1]]) %>% 
  bind_rows(data.frame(eBird2016$data[[1]])) %>% 
  bind_rows(data.frame(eBird2017$data[[1]]))


eBird <- st_as_sf(x = eBird,
                  coords = c('longitude', 'latitude'),
                  crs = proj)

eBird$Year <- eBird$year

eBird <- eBird[unlist(st_intersects(AL, eBird)),]

We then get onto the model description, which in this case includes a shared spatial field between the two datasets. This shared spatial field is characterized by an ar1 process. To add this structure into the model, we specify the parameter temporalModel in the function intModel appropriately, Furthermore we specified the hyper parameters for both the random field and the temporal effect.


hyperParams <- list(rho = list(prior = "pc.prec", param = c(0.01, 0.01)))

modelSetup <- intModel(eBird, BBSColinusVirginianus,
                       Coordinates = c('Longitude', 'Latitude'),                                                    temporalName = 'Year',
                       Projection =  proj, Mesh = mesh, 
                       responsePA =  'NPres', trialsPA = 'Ntrials',
                       temporalModel = list(model = 'ar1', hyper = hyperParams))

modelSetup$specifySpatial(sharedSpatial = TRUE, prior.sigma = c(0.2, 0.01), 
                          prior.range = c(0.4, 0.01))


modelSetup$plot()


modelSetup$changeComponents()

Next we run the model, using the function fitISDM. Due to time considerations, inference for this model is not completed in the vignette. However, both the data and script is provided for the user to complete the analysis.


mod <- fitISDM(modelSetup,
                options = list(control.inla = list(int.strategy = 'eb')))


preds <- predict(mod, mask = AL, mesh = mesh, temporal = TRUE, fun = '')

plot_preds <- plot(preds, whattoplot = 'median', plot = FALSE)

plot_preds + 
  geom_sf(data = st_boundary(AL), lwd = 1.2) + 
  scico::scale_fill_scico(palette = "lajolla") + 
  theme_minimal()

Spatiotemporal example

Philip Mostert

2024-02-02