Experiments with sdmTMB models • sdmTMBexperiments

sdmTMBexperiments collects experimental helpers for fitted sdmTMB models. This vignette walks through the main functions, fitting small/fast models throughout so the whole vignette runs.

library(sdmTMBexperiments)
library(sdmTMB)
library(ggOceanMaps)

Model comparison

Fit a couple of models and compare their indices and spatial predictions. Mohn’s rho indicates the summed proportional difference during the final year and makes only sense when comparing multiple models. Included here as an example. Note that the plot() functions return standard ggplot2 figures and you can add layers or modify them as you would any other ggplot object.

mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 20)
m1 <- sdmTMB(
  density ~ 0 + as.factor(year),
  data = pcod,
  time = "year",
  mesh = mesh,
  family = tweedie()
)
m2 <- sdmTMB(
  density ~ 0 + as.factor(year) + depth_scaled,
  data = pcod,
  time = "year",
  mesh = mesh,
  family = tweedie()
)
nd <- replicate_df(qcs_grid, "year", unique(pcod$year))

cmp <- compare_models(m1, m2, newdata = nd, object_crs = 32609)
#> ℹ `ln_tau_O` is an internal parameter affecting `sigma_O`
#> ℹ `sigma_O` is the spatial standard deviation
#> ℹ `ln_tau_E` is an internal parameter affecting `sigma_E`
#> ℹ `sigma_E` is the spatiotemporal standard deviation
#> ℹ `ln_kappa` is an internal parameter affecting `range`
#> ℹ `range` is the distance at which data are effectively independent

# Index time series, annotated with Mohn's rho and modified using ggplot2 syntax
plot(cmp, mohn = TRUE) +
  ggplot2::scale_x_continuous(
    breaks = seq(min(pcod$year), max(pcod$year), by = 2),
    labels = substr(seq(min(pcod$year), max(pcod$year), by = 2), 3, 4)
  )
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.


# per-cell prediction-difference maps
plot(cmp, type = "spatial", fill_prefix = "Biomass density")

The spatial_response argument controls which per-cell summary of the prediction difference is mapped, since a model comparison typically spans several years and these need to be summarised down to one map. All three options are computed from the response-scale (not log-scale) predictions. "mean" is the average absolute difference across years, in the response’s own units (e.g. density); it is not comparable between cells with very different baseline abundance. "mean_pr_diff" (the default, used above) is the percent difference between the across-year mean predictions — a ratio of averages, not an average of ratios. "pr_mean" is instead the average of each year’s own percent difference; because each year is divided by that year’s own reference value, a handful of low-abundance years/cells can dominate this average with large, noisy percentages, which is why "mean_pr_diff" is the more robust default. See ?plot.sdmTMB_model_comparison for the exact formulas.

In practice, "pr_mean" and "mean_pr_diff" are usually close to each other and can look nearly identical on the map: m1 and m2 differ by a single covariate here, and the reference model’s predictions don’t swing wildly from year to year within a cell, so a mean-of-ratios and a ratio-of-means end up similar. They are not identical, though — the two summaries are highly correlated but diverge most for cells that have at least one year with an unusually small reference-model estimate, which inflates that single year’s percent difference in "pr_mean" (but not "mean_pr_diff", which only sees the multi-year average).

# Aggregate to one value per cell, exactly as plot(..., type = "spatial") does
agg <- cmp$grid[[1]] |>
  dplyr::group_by(X, Y) |>
  dplyr::reframe(
    pr_mean = mean(pr_diff, na.rm = TRUE),
    mean_pr_diff = mean(mean_pr_diff, na.rm = TRUE)
  )
cor(agg$pr_mean, agg$mean_pr_diff)
#> [1] 0.9993678
max(abs(agg$pr_mean - agg$mean_pr_diff))
#> [1] 57.80477

plot(
  cmp,
  type = "spatial",
  fill_prefix = "Biomass density",
  spatial_response = "pr_mean"
)

Diagnostics and result maps

All map helpers render on a ggOceanMaps::basemap() when ggOceanMaps is installed (otherwise on a plain map). Plotting the SPDE mesh and model results requires knowing the utm_crs used by the sdmTMB::add_utm_columns() function. The X and Y coordinates are divided by 1000 in that function. The example dataset does not contain information about the UTM zone, but the information is extracted from sdmTMB examples. In real data, use the sdmTMB::get_crs() function to get the correct UTM zone.

plot_mesh(mesh, crs = 32609, data = pcod, value = "density")


pred <- predict(m1, newdata = nd, type = "response") # response scale, not link scale
plot_prediction_map(pred, crs = 32609, transform = "sqrt") # density averaged over years

plot_random_field(pred, crs = 32609, field = "omega") # spatial field

plot_residual_qq(m1)

plot_variance_partitioning(m1)

Model diagnostics plots automatically detect delta models and plot all components. The example below fits a delta model and plots the residuals and variance partitioning.

m3 <- sdmTMB(
  density ~ 0 + as.factor(year),
  data = pcod,
  time = "year",
  mesh = mesh,
  family = delta_lognormal()
)

plot_residual_qq(m3)

plot_variance_partitioning(m3)

Mesh-resolution sensitivity

Mesh resolution can affect model results. compare_meshes() fits a model at multiple mesh resolutions and compares the resulting indices and spatial predictions. The example below compares two meshes with cutoffs of 10 and 30 km, using the 20 km mesh as the reference, and shows how finer mesh size leads to generally lower index estimates.

mc <- compare_meshes(
  m2,
  mesh_cutoffs = c(10, 30),
  newdata = nd,
  object_crs = 32609,
  reference_name = "20"
)
#> Model for cutoff 30 did not pass sdmTMB::sanity() checks.
#> Model for cutoff 20 did not pass sdmTMB::sanity() checks.
plot(mc, cutoff_label = "Cutoff (km)")

plot(mc, type = "spatial", spatial_response = "mean")

The sdmTMB_mesh_cv() performs a block cross-validation of a fitted model at multiple mesh resolutions. The example below compares the cutoff values in the example above and indicates that the coarsest mesh would provide the best model fit, which is somewhat surprising. The function is still very experimental, however. The blockCV package is required for this function.

# Requires the 'blockCV' package.
cv <- sdmTMB_mesh_cv(m2, mesh_cutoffs = c(10, 20, 30), k = 3, object_crs = 32609)
cv
#> # A tibble: 3 × 10
#>   cutoff crashed converged sanity sum_loglik mean_aic deltaNLL deltaNLL_pr
#>    <dbl> <lgl>   <lgl>     <lgl>       <dbl>    <dbl>    <dbl>       <dbl>
#> 1     10 FALSE   TRUE      TRUE       -6705.    8561.       0         0   
#> 2     20 FALSE   TRUE      FALSE      -6845.    8676.     140.        2.08
#> 3     30 FALSE   TRUE      FALSE      -7089.    8752.     383.        5.72
#> # ℹ 2 more variables: deltaAIC <dbl>, deltaAIC_pr <dbl>
plot(cv)

Retrospective analysis

Retrospective analysis may be required in assessment context to demonstrate model stability. The sdmTMB_retro() function fits a model to the full time series, then sequentially removes the most recent years and refits the model. The example below removes the last 5 years of data and plots the resulting indices. The index_divisor and index_unit arguments are used to scale the index for plotting.

retr <- sdmTMB_retro(m1, nd, nyears = 5)
plot(retr, index_divisor = 1e3, index_unit = "1000 t")

Prediction grids

The bundled sebastes dataset (beaked redfish survey stations) is used to demonstrate the prediction-grid helpers. make_survey_prediction_region() builds a concave hull around a survey’s stations, and make_prediction_grid() lays a regular grid over it. The function can be used to extract depths and drop cells outside a specified depth range. The example below builds a prediction region and grid for the EggaN survey series, and extracts depths. By default depths come from ggOceanMaps::get_depth() (which downloads the ETOPO 2022 bathymetry on first use); here we pass depth_fun = sebastes_depth, a small bundled bathymetry of the study area, so the vignette needs no download. Use the default for real analyses.

data(sebastes)
eggan <- subset(sebastes, cruiseseries == "EggaN")

region <- make_survey_prediction_region(eggan, crs = 32633, concavity = 5)

# the result is an sf object
qmap(region, fill = I(NA))


# max_depth turns on depth extraction, via ggOceanMaps::get_depth() by
# default (bathy_style controls the bathymetry source); cells outside
# c(1, max_depth) are dropped. plot() then colours by depth automatically.
# Here `depth_fun = sebastes_depth` uses a small bundled demonstration
# bathymetry to avoid the large download get_depth() needs on first use.
grid <- make_prediction_grid(
  region, resolution = 20, max_depth = 1000, depth_fun = sebastes_depth
)
range(grid$depth)
#> [1]  10 970
plot(grid)

Survey comparison and prediction grids

When compiling a survey index from multiple survey time-series, it is often useful to compare the survey-specific indices to the compiled index to understand which surveys carry the signal in the compiled index. The compare_surveys() function fits a model per survey series, builds each survey’s prediction grid, and compares the survey indices to a compiled reference index. The example below uses the bundled sebastes dataset (beaked redfish survey stations) to demonstrate this functionality. EcoS and WinterS are the two series with a coverage over the Barents Sea, so we use those two, restricted to a few recent years and a coarse mesh to keep this vignette fast.

Both compiled_model and the survey models below are fit with a swept_area offset, so the response is biomass (raw catch), not density: density is already swept_area-normalised (density = biomass / swept_area), so fitting it together with a swept_area offset would double-correct for effort. Predicting onto full_grid (which has no swept-area column) makes predict.sdmTMB() print an informational message that it is assuming an offset of zero for the prediction — this is expected: grid cells aren’t real trawl hauls, so we deliberately predict at a zero offset (i.e. per 1 swept-area unit, which is exactly density’s definition) and let the area argument of get_index_split() / compare_surveys() do the area weighting instead. compare_surveys() / compare_models() / compare_meshes() already suppress this message internally, so it won’t show when calling those; it can still show for a direct get_index_split()/predict() call on an offset model, such as compiled_index below. swept_area is recorded in square nautical miles (see ?sebastes), not the km² implied by the X/Y grid coordinates, so cell areas must be converted to nautical miles before squaring.

sub <- subset(sebastes, cruiseseries %in% c("EcoS", "WinterS") & year >= 2020)

# A simple "compiled" full-domain model and index, used as the reference
compiled_mesh <- make_mesh(sub, c("X", "Y"), cutoff = 50)
compiled_model <- sdmTMB(
  biomass ~ 0 + as.factor(year),
  data = sub,
  time = "year",
  mesh = compiled_mesh,
  family = tweedie(),
  spatiotemporal = "iid",
  offset = log(sub$swept_area)
)

region <- make_survey_prediction_region(sub, crs = 32633, concavity = 5)
grid <- make_prediction_grid(
  region, resolution = 40, max_depth = 500, depth_fun = sebastes_depth
)

plot(grid)


full_grid <- replicate_df(
  grid,
  "year",
  unique(sub$year)
)

# swept_area is in square nautical miles: convert the 40 km grid resolution
# to nautical miles before squaring, so cell area matches those units.
cell_area <- (40 / 1.852)^2

compiled_index <- get_index_split(
  compiled_model,
  newdata = full_grid,
  area = cell_area,
  silent = TRUE
)

sc <- compare_surveys(
  data = sub,
  compiled_index = compiled_index,
  compiled_model = compiled_model,
  surveys = c("EcoS", "WinterS"),
  newdata = full_grid,
  object_crs = 32633,
  mesh_cutoff = 50,
  offset_col = "swept_area",
  area = cell_area,
  verbose = FALSE
)

plot(sc, index_divisor = 1e6, unit = "kt") # scaled index time series

plot(sc, type = "scatter") # survey vs compiled

The two prediction-difference maps are tall relative to their width (the survey domain’s north-south extent is larger than its east-west extent), so this last plot uses a taller fig.height than the defaults above to fit both side by side without wasted space.

plot(sc, type = "spatial", spatial_response = "mean") # prediction differences

The individual survey indices are lower than the compiled index because the survey models are fit to a smaller spatial domain than the compiled model. ### Survey-specific settings and the same-area comparison

Surveys differ in coverage, season, depth range and sampling frequency, so a model that suits the compiled data set rarely suits each survey on its own. survey_spec() gives a survey its own series, model settings and prediction grid. A survey can combine several series: two series that alternate years, for example, can be fitted together with the series as a fixed effect, AR(1) spatiotemporal fields that bridge the gaps, and predictions for one of the series (predict_values). depth_range = "data" limits the prediction grid to the 2–98% quantiles of the survey’s haul depths. Settings not given in survey_spec() are taken from the compare_surveys() arguments, which in turn are inherited from the compiled model.

Because each survey index covers only its own area, it is not directly comparable with the compiled index. With same_area = TRUE (the default when compiled_model is given), both the survey model and the compiled model are also predicted over the part of the survey grid inside the compiled domain. plot(type = "same_area") then shows two estimates of the same quantity for the same area: where they differ, the compiled model, informed by all surveys, departs from what the survey alone says about its own area. sc$settings lists the settings used for each survey.

sc2 <- compare_surveys(
  data = sub,
  compiled_index = compiled_index,
  compiled_model = compiled_model,
  surveys = list(
    EcoS = "EcoS",
    "EcoS + WinterS" = survey_spec(
      c("EcoS", "WinterS"),
      formula = biomass ~ 0 + cruiseseries,
      spatiotemporal = "ar1",
      predict_values = list(cruiseseries = "EcoS"),
      depth_range = "data"
    )
  ),
  newdata = full_grid,
  object_crs = 32633,
  mesh_cutoff = 50,
  offset_col = "swept_area",
  area = cell_area,
  verbose = FALSE
)

sc2$settings[c("survey", "series", "spatiotemporal", "min_depth", "max_depth", "grid_cells")]
#> # A tibble: 2 × 6
#>   survey         series        spatiotemporal min_depth max_depth grid_cells
#>   <chr>          <chr>         <chr>              <dbl>     <dbl>      <int>
#> 1 EcoS           EcoS          iid                 NA         NA         400
#> 2 EcoS + WinterS EcoS, WinterS ar1                 63.2      573.        413
plot(sc2, type = "same_area", index_divisor = 1e6, unit = "kt")

Which survey carries the signal?

sdmTMB_leave_out() refits a model once without each survey series (or cruise, or haul) and compares the resulting index with a refit to all data. Years with no data left after removing a series are kept in the refit as extra_time and reported as NA. A large change when a series is removed shows that the series drives the compiled index in those years. index = "split" returns confidence intervals, which suits a few identifiers such as survey series; the default plug-in index is faster for screening many hauls.

lo <- sdmTMB_leave_out(
  compiled_model,
  id_col = "cruiseseries",
  newdata = full_grid,
  area = cell_area,
  index = "split",
  verbose = FALSE
)

lo$influence
#> # A tibble: 2 × 9
#>   id      n_removed id_years  haul_share influence leverage max_abs_change
#>   <chr>       <int> <list>         <dbl>     <dbl>    <dbl>          <dbl>
#> 1 WinterS      2016 <int [6]>      0.656     0.138    0.211          0.277
#> 2 EcoS         1056 <int [6]>      0.344    -0.117   -0.341          0.197
#> # ℹ 2 more variables: year_max_change <int>, all_ok <lgl>
plot(lo, index_divisor = 1e6, unit = "kt")