In pbs-assess/sdmTMB: Spatial and Spatiotemporal SPDE-Based GLMMs with 'TMB'

If the code in this vignette has not been evaluated, a rendered version is available on the documentation site under 'Articles'.

dplyr_installed <- require("dplyr", quietly = TRUE)
ggplot_installed <- require("ggplot2", quietly = TRUE)
pkgs <- dplyr_installed && ggplot_installed
EVAL <- identical(Sys.getenv("NOT_CRAN"), "true") && pkgs
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.asp = 0.618,
  eval = EVAL,
  purl = EVAL
)

library(ggplot2)
library(dplyr)
library(sdmTMB)

Using the built-in British Columbia Queen Charlotte Sound Pacific Cod dataset, we might be interested in fitting a model that describes spatially varying trends through time. The data are as follows:

• There are columns for depth and depth squared.
• Depth was centred and scaled by its standard deviation and we've included those in the data frame so that they could be used to similarly scale the prediction grid.
• The density units should be kg/km^2^.
• Here, X and Y are coordinates in UTM zone 9.

pcod_spde <- make_mesh(pcod, c("X", "Y"), cutoff = 5)
plot(pcod_spde)

We will fit a model that includes an overall intercept, a slope for 'year', an intercept spatial random field, and another random field for spatially varying slopes that represent temporal trends that vary spatially (spatial_varying argument).

First, we will set up a column for time that is centered to help break correlation with the intercept and help estimation (and interpretation of the intercept) and we will scale the year values so that a unit of change is a decade instead of a year. This is likely more interpretable and can help with estimation by keeping the coefficients on a reasonable scale.

d <- pcod
d$scaled_year <- (pcod$year - mean(pcod$year)) / 10
unique(d$scaled_year)

Now we will fit a model using spatial_varying ~ 0 + scaled_year. The 0 + drops the intercept, which is already present due to spatial = "on", which is the default.

fit <- sdmTMB(
  density ~ scaled_year, 
  data = d,
  mesh = pcod_spde, 
  family = tweedie(link = "log"),
  spatial = "on",
  spatial_varying = ~ 0 + scaled_year
)

print(fit)
sanity(fit)

We have not included spatiotemporal random fields for this example for simplicity, but they could also be included.

Let's extract some parameter estimates. Look for sigma_Z coefficients:

tidy(fit, conf.int = TRUE)
tidy(fit, "ran_pars", conf.int = TRUE)

Let's look at the predictions and estimates of the spatially varying coefficients on a grid. First, we will create a small function to help with plotting:

plot_map_raster <- function(dat, column = est) {
  ggplot(dat, aes(X, Y, fill = {{ column }})) +
    geom_raster() +
    coord_fixed() +
    scale_fill_gradient2()
}

We need to predict on a grid to make our visualizations. We also need to add a column for scaled_year to match the fitting. Make sure you scale based on the same values! Here, that means using the mean from our fitted data and turning years into decades.

nd <- replicate_df(qcs_grid, "year", unique(pcod$year))
nd$scaled_year <- (nd$year - mean(pcod$year)) / 10
pred <- predict(fit, newdata = nd)

Let's look at the spatial trends. The zeta_s_scaled_year in the predictions are values from a random field describing deviations from the overall effect of year_scaled on species density. These are in link (log) space. These are the spatially varying coefficients.

plot_map_raster(pred, zeta_s_scaled_year)

If we wanted to get at the actual slope of scaled_year in space, we would also have to add on the main overall (fixed effect) coefficient.

coefs <- tidy(fit, conf.int = TRUE)
scaled_year_coef <- coefs$estimate[coefs$term == "scaled_year"]
scaled_year_coef
exp(scaled_year_coef)

Our main effect tells us that the overall linear trend of density has been a multiplicative factor of r round(exp(scaled_year_coef), 2) per decade. Conversely, the overall trend has been a decline of r 100 * round(1 - exp(scaled_year_coef), 2)% per decade.

We can add this overall effect to the spatial deviations from this effect to get the spatially varying trend:

plot_map_raster(pred, scaled_year_coef + zeta_s_scaled_year)

We could put those slopes back into natural space:

plot_map_raster(pred, exp(scaled_year_coef + zeta_s_scaled_year)) +
  scale_fill_gradient2(midpoint = 1)

Here, 2, for example, means a 2-fold change in density in that location per decade.

We might prefer to log transform the color scale. This gets us back to the plot in log space but with values on the colour legend that represent (multiplicative) values in natural space.

plot_map_raster(pred, exp(scaled_year_coef + zeta_s_scaled_year)) +
  scale_fill_gradient2(midpoint = 0, trans = "log10")

The fastest way to get uncertainty on this spatially varying coefficient is with draws from the joint precision matrix of the parameters:

set.seed(82938)
psim <- predict(fit, newdata = nd, nsim = 200, sims_var = "zeta_s")
pred$se <- apply(psim, 1, sd)
plot_map_raster(pred, se) +
  scale_fill_viridis_c()

We can see a visual pattern of the SPDE mesh in this plot. The uncertainty standard error is most accurate at the nodes (vertices, knots) of the mesh.

pbs-assess/sdmTMB documentation built on May 17, 2024, 11:31 a.m.