INLAutils (version 0.0.5)

inlaSDM: Fit INLA species distribution models.

Description

Fit INLA species distribution models.

Usage

inlaSDM(dataframe, predictors, include = 1:raster::nlayers(predictors),
  step = FALSE, invariant = "0 + Intercept", cross_validation = FALSE,
  cv_folds = 5, spatial = TRUE, num.threads = 1,
  meshvals = list(inner.max.edge = max(raster::res(predictors)) * 10,
  outer.max.edge = max(raster::res(predictors)) * 100, cutoff = 0, inner.offset
  = -0.1, outer.offset = -0.3))

Arguments

dataframe

A SpatialPointsDataFrame containing the presence absence values

predictors

Raster of predictors (covariates)

include

Vector of integers describing which covariates to include in the model

step

Logical indicating whether to run stepwise elimination of variables.

invariant

Character indicating the parts of the model formula that should not change despite stepwise selection (e.g. the intercept).

cross_validation

Run cross validation?

cv_folds

How many folds should the data be split into?

spatial

Run INLA with a spatial term.

num.threads

How many threads should be used for parallel computation.

meshvals

List giving details for the creation of the INLA mesh (see details and inla.mesh.2d)

Details

For now invariant MUST include 'Intercept'.

meshvals takes a list of up to five named values:

inner.max.edge

Maximum triangle length for inner domain.

outer.max.edge

Maximum triangle length for outer domain.

cutoff

Minumum allowed distance between mesh nodes.

inner.offset

Extension distance beyond points.

outer.offset

Additional extension distance with larger triangles (with max length outer.max.edge)

Note that negative values for the offsets are in absolute units by default. Negative values give the extension distance relative to the diameter of the coordinates range (i.e. -0.1 will create an extension 10% the that 10% the diameter of the points).

These values are explained in more detail in inla.mesh.2d.

Examples

# NOT RUN {
library(INLA)
set.seed(6)

# Create locations, presence absence points and covariates 
#   with spatial and environmental relationships
coords <- data.frame(long = c(rnorm(70), rnorm(30, 3)), lat = rnorm(100))
PA <- rep(c(0, 1), each = 50)
x <- data.frame(x1 = rnorm(100), # no relationship
               x2 = c(rnorm(70), rnorm(30, 5))) # positive relationship

# Have a look
\dontrun{
ggplot(cbind(x, PA), aes(x1, PA)) + 
 geom_point() +
 geom_smooth(method = 'glm', method.args = list(family = 'binomial'))


ggplot(cbind(x, PA), aes(x2, PA)) + 
 geom_point() +
 geom_smooth(method = 'glm', method.args = list(family = 'binomial'))


ggplot(cbind(coords, PA), aes(long, lat, colour = PA)) + geom_point()
}

# Set raster resolution
res <- 50

# Create raster limits
xrange <- range(coords$long)
xrange <- c(floor(xrange[1]), ceiling(xrange[2]))
yrange <- range(coords$lat)
yrange <- c(floor(yrange[1]), ceiling(yrange[2]))

# Calculate number of cells
xcells <- res * (xrange[2] - xrange[1])
ycells <- res * (yrange[2] - yrange[1])

# Create an empty raster of correct dims
suppressWarnings(
 raster <- raster::raster(matrix(NA, ncol = ycells, nrow = xcells), 
   xmn = xrange[1], xmx = xrange[2], ymn = yrange[1], ymx = yrange[2])
)
# Add dataframe data to rasters, then fill gaps with random data.
x1 <- raster::rasterize(coords, raster, x$x1)
x1[is.na(x1)] <- rnorm(sum(is.na(raster::getValues(x1))))

x2 <- raster::rasterize(coords, raster, x$x2)
x2[is.na(x2)] <- rnorm(sum(is.na(raster::getValues(x2))))

# Stack rasters
predictors <- raster::stack(x1, x2)

# Pull together coordinates and PA data into SpatialPointsDataFrame
dataframe = sp::SpatialPointsDataFrame(coords = coords, data = data.frame(y = PA))

# Run the model.
model <- inlaSDM(dataframe, 
                 predictors, 
                 spatial = TRUE, 
                 cross_validation = FALSE,
                 meshvals = list(cutoff = 0.3, inner.max.edge = 1))
                 
autoplot(model$mesh[[1]])
autoplot(model$models[[1]])
# }