inlaSDM: Fit INLA species distribution models.
In INLAutils: Utility Functions for 'INLA'

Description Usage Arguments Details Examples

Fit INLA species distribution models.

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))

`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`)

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.

## 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]])

## End(Not run)