Select area

Next we specify the study domain for the study: in this case Norway. This can be achieved using either the countryName argument which will then access the object from the giscoR (Hernangómez 2023) R package, or by using Object and supplying our own polygon object. In this case we choose to add our own polygon object because the rigid coastline of a unedited Norway polygon will cause issues further down.

We then use .$plot() to see what the boundary looks like.

yes <- FALSE
if (yes) {
Norway <- giscoR::gisco_get_countries(country = 'Norway', resolution = 60)
Norway <- st_cast(st_as_sf(Norway), 'POLYGON')
Norway <- Norway[which.max(st_area(Norway)),]
Norway <- rmapshaper::ms_simplify(Norway, keep = 0.8)

workflow$addArea(Object = Norway)
workflow$plot()
}

Norway <- readRDS('IntegratedLakefish/Norway.rds')
workflow$addArea(Object = Norway)

Adding occurrence data

Species’ occurrence data is certainly the most important component of a SDM, and intSDM has two slot functions to help you add data into the workflow: .$addGBIF() and .$addStructured().

The former of which uses the rgbif package to download data directly from GBIF. For this function, we need to specify the name of the dataset (datasetName) and the type of the dataset (datasetType) – which can be one of PO, PA or Counts. The … argument is used to specify any addiditional arguments for rgbif::occ_data() (Chamberlain, Oldoni, and Waller 2022) (in this case, limit and datasetKey). If datasetType = 'PA', absences may be generated using generateAbsences = TRUE. This will treat the obtained data as a checklist survey data: combining all the sampling locations for the species in the dataset, and creating absences when a given species did not occur in a given region.

For this example we consider three sources of data. The Vascular Plant Field Notes is a collection of observations provided by the Norwegian University of Science and Technology’s (NTNU) (Norwegian University of Science and Technology 2023) University Museum and the University of Oslo (UiO) (University of Oslo 2023), containing records of standardized cross-lists of most vascular plants found in Norway. We treat these two datasets as detection/non-detection data, generating absences in sampling locations where the species does not occur.

The other source of data considered comes from the Norwegian Species Observation service (published by Artsdatabanken) (The Norwegian Biodiversity Information Centre 2023). This data is a collection of citizen science records – and as a result we treat it as presence-only data.

workflow$addGBIF(datasetName = 'NTNU', 
                 datasetType = 'PA',
                 limit = 10000,
                 coordinateUncertaintyInMeters = '0,50',
                 generateAbsences = TRUE, 
                 datasetKey = 'd29d79fd-2dc4-4ef5-89b8-cdf66994de0d')

workflow$addGBIF(datasetName = 'UiO',
                 datasetType = 'PA',
                 limit = 10000,
                 coordinateUncertaintyInMeters = '0,50',
                 generateAbsences = TRUE,
                 datasetKey = 'e45c7d91-81c6-4455-86e3-2965a5739b1f')

workflow$addGBIF(datasetName = 'CZ', 
                 datasetType = 'PO',
                 coordinateUncertaintyInMeters = '0,50',
                 limit = 10000,
                 datasetKey = 'b124e1e0-4755-430f-9eab-894f25a9b59c')

workflow$plot(Species = TRUE)

Adding covariate data

Covariate data may be added to the model using .$addCovariates(). Layers from WorldClim (Fick and Hijmans 2017) may be accessed using the worldClim argument. This in turn uses the geodata (Hijmans et al. 2023) R package to obtain spatRaster objects of the covariates cropped around the study domain. Other covariate layers may be added using the Object argument.

workflow$addCovariates(worldClim = 'tavg', res = 5, Function = scale)
workflow$plot(Covariates = TRUE)

Metadata

We can then view the metadata for the obtained occurrence records using the .$obtainMeta() function, which will give us the citation for the datasets used in this workflow.

workflow$obtainMeta()

Creating an inla.mesh object

One of the objects required for our model is an inla.mesh object, which we will use in the approximation of our spatial random fields. The .$addMesh() function’s argument … uses INLA::inla.mesh.2d() to create this object.

workflow$addMesh(cutoff = 20 * 0.25,
                 max.edge = c(60, 80)*0.5, #0.25
                 offset= c(30, 40))

workflow$plot(Mesh = TRUE)

Specify priors

Furthermore we also used penalizing complexity (PC) priors in our model, which are designed to control the spatial range and standard deviation in the GRF’s Matérn covariance function in order to reduce over-fitting in the model (Simpson et al. 2017).

workflow$specifySpatial(prior.range = c(100, 0.1),
                        prior.sigma = c(1, 0.1),
                        constr = FALSE)

We also specify priors for the intercept terms and the fixed effects of the models. In this vignette we choose tight priors (mean = 0; precision = 1).

workflow$specifyPriors(effectNames = 'Intercept',
                       Mean = 0, Precision = 1)

workflow$specifyPriors('tavg', Mean = 0, Precision = 1)

Bias field

We specify an additional spatial effect for the citizen science data using .$biasFields() to account biases in the collection process (Simmonds et al. 2020).

workflow$biasFields('CZ', 
                    prior.range = c(100, 0.1), #1 #0.1
                    prior.sigma = c(1, 0.1), #1 #0.1
                    constr = FALSE)

Model options

We specify the model output using the function .$workflowOutput. In this workflow, we want to return the R-INLA model objects and the maps of the predicted intensity.

#workflow$crossValidation(Method = 'Loo')
workflow$workflowOutput(c('Maps', 'Model', 'Bias'))