Compare individual survey indices to a compiled sdmTMB index

Fits separate sdmTMB models for individual surveys and compares the resulting indices to a compiled reference index. A survey can be one survey series or several series combined, and each survey can have its own model settings and prediction grid (see survey_spec()).

Two comparisons are made:

Survey area: each survey index is predicted over the survey's own footprint (a concave hull of its hauls, make_survey_prediction_region(), optionally limited to a depth band) and plotted against the compiled index for the whole domain (plot(type = "index")). The two cover different areas, so only their trends are comparable.
Same area (same_area = TRUE): the survey model and the compiled model are both predicted over the part of the survey grid inside the compiled prediction domain, giving two indices of the same quantity for the same area (plot(type = "same_area")). Differences show where the compiled model, informed by all surveys, departs from what the survey alone says about its own area.

To find which survey carries the signal in the compiled index, refit the compiled model without each survey with sdmTMB_leave_out().

Usage

compare_surveys(
  data,
  compiled_index,
  compiled_model = NULL,
  surveys = NULL,
  newdata,
  object_crs = NULL,
  formula = NULL,
  family = NULL,
  spatial = NULL,
  spatiotemporal = NULL,
  control = NULL,
  mesh_cutoff = 30,
  offset_col = NULL,
  survey_col = "cruiseseries",
  xy_cols = c("X", "Y"),
  coord_multiplier = 1000,
  bias_correct = FALSE,
  level = 0.95,
  area = 1,
  silent = TRUE,
  nsplit = NULL,
  n_workers = 1,
  future_seed = TRUE,
  concavity = 5,
  buffer = 1,
  depth_range = NULL,
  depth_quantiles = c(0.02, 0.98),
  depth_col = "depth",
  same_area = !is.null(compiled_model),
  compiled_region = NULL,
  keep_models = FALSE,
  reference_label = "Compiled",
  verbose = TRUE
)

Arguments

data

A data frame containing the modelling data, with a column identifying the survey series (see survey_col) and coordinate columns (see xy_cols).

compiled_index

A data frame of the compiled reference index, with columns year, est, lwr, upr, se. A time column named after compiled_model$time (as returned by sdmTMB::get_index_split()) is renamed to year. It must be in the same units as the survey indices, that is, sdmTMB::get_index_split() output for newdata and area without rescaling; a warning is given when the two differ by more than 1000 times.

compiled_model

A fitted sdmTMB::sdmTMB() model for the compiled index, or NULL (default). When provided, all settings of its original call (as with sdmTMB::update.sdmTMB(); e.g. family, random fields, time column, spatial_varying, priors, ...), its coordinate columns and CRS are inherited, except those overridden by formula, family, spatial, spatiotemporal, control and offset_col, or by a survey_spec(). Its data and mesh are replaced by each survey's data and a new mesh. A numeric offset or weights vector in the original call cannot be carried over to the survey subsets and gives an error (use offset_col or column-name arguments instead). The same-area comparison and spatial prediction differences are computed. When NULL, these are skipped, formula, family and object_crs must be supplied, and the survey models use only formula, family, spatial, spatiotemporal, control, offset_col, a "year" time column and the xy_cols coordinates (plus any survey_spec() settings).

surveys

Surveys to compare. One of:

NULL (default): every unique value of data[[survey_col]];
a character vector of series names (names, if given, are used as labels);
a named list whose elements are character vectors of series combined into one survey index, or survey_spec() objects with survey-specific settings, e.g. list(EcoS = "EcoS", "Egga N+S" = survey_spec(c("EggaN", "EggaS"), spatiotemporal = "ar1")).

newdata

A data frame (the compiled prediction grid, with a time column) from which survey grids are clipped. Also defines the compiled domain for the same-area comparison.

object_crs

Coordinate reference system. If NULL, taken from compiled_model.

formula

Model formula. If NULL, inherited from compiled_model.

family

sdmTMB family. If NULL, inherited from compiled_model.

spatial, spatiotemporal

Random field settings. If NULL, inherited from compiled_model (falling back to "on"/"ar1").

control

An sdmTMB::sdmTMBcontrol() object. If NULL, inherited from compiled_model (falling back to sdmTMBcontrol(newton_loops = 2)).

mesh_cutoff

Numeric. Default mesh cutoff (in xy_cols units, typically km) for the survey models. Default 30.

offset_col

Character. Column in data to use as a log-offset (e.g. "swept_area"). NULL (default) means no offset.

survey_col

Character. Column identifying the survey series. Default "cruiseseries".

xy_cols

Character vector of the coordinate columns in data and newdata. Only used when compiled_model = NULL; otherwise taken from the model mesh. Default c("X", "Y").

coord_multiplier

Numeric factor to convert xy_cols to CRS units, used to build and clip to the survey footprints and to draw the spatial maps. Default 1000 (coordinates in km, CRS in m); use 1 for coordinates already in CRS units.

bias_correct

Logical. Apply bias correction? Default FALSE.

level

Numeric. Confidence level for index intervals. Default 0.95.

area

Grid cell area for index calculation: a single value, a column name in newdata, or a numeric vector with one value per row of newdata (subset to each survey's clipped grid). Passed to sdmTMB::get_index_split().

silent

Logical. Suppress model output? Default TRUE.

nsplit

Integer. Number of splits for sdmTMB::get_index_split(). NULL (default) uses about 75 000 grid rows per split, which avoids the silent crashes of bias-corrected extraction on very large grids.

n_workers

Integer. Number of parallel workers (surveys are processed in parallel). Default 1, which uses whatever future::plan() is currently set. Values > 1 set a temporary parallel plan (forked future::multicore where supported, otherwise future::multisession).

future_seed

Logical or integer. Future-compatible random seed.

concavity

Numeric. Default concavity for the survey footprint hulls (see make_survey_prediction_region()). Default 5.

buffer

Numeric. Default buffer around the survey footprint hulls, in CRS units (typically m). Default 1.

depth_range

Default depth band for the survey grids: NULL (no depth limit, the default), a numeric vector c(min, max), or "data" for quantiles of each survey's haul depths.

depth_quantiles

Numeric vector of two probabilities used when depth_range = "data". Default c(0.02, 0.98).

depth_col

Character. Depth column in data and the prediction grids. Default "depth".

same_area

Logical. Compute the same-area comparison? Default TRUE when compiled_model is given.

compiled_region

Optional sf polygon of the compiled prediction domain, used to limit custom survey grids (survey_spec() newdata) to the compiled domain for the same-area comparison. If NULL (default), a concave hull of the unique newdata coordinates is used. Grids clipped from newdata are inside the domain by construction.

keep_models

Logical. Return the fitted survey models? Default FALSE; models can be large.

reference_label

Character. Label for the compiled reference index. Default "Compiled".

verbose

Logical. Show progress messages (via cli when installed)? Default TRUE.

Value

An object of class sdmTMB_survey_comparison, a list with

sanity: sdmTMB::sanity() results, one column per survey;
settings: one row per survey with the series, number of hauls, years, depth band, mesh size, number of grid cells and model settings;
index: the indices, with columns model (survey label or reference_label), survey, domain ("full" for the compiled index over the whole domain, "survey" for a survey index over its own footprint, "shared" for both models over the same area), year, est, lwr, upr, se, cv. The time column is always named year;
grid: per-cell prediction differences between each survey model and the compiled model over the shared area (or NULL);
models: the fitted survey models if keep_models = TRUE.

Carries crs, coord_multiplier, reference_label and surveys attributes.

Author

Mikko Vihtakari

Examples

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

# A simple "compiled" full-domain model and index, used as the reference.
# `biomass` (raw catch), not `density`, is the response here: `density` is
# already swept_area-normalised, so using it together with the swept_area
# offset below would double-correct for effort.
compiled_mesh <- sdmTMB::make_mesh(sub, c("X", "Y"), cutoff = 50)
compiled_model <- sdmTMB::sdmTMB(
  biomass ~ 0 + as.factor(year), data = sub, time = "year",
  mesh = compiled_mesh, family = sdmTMB::tweedie(), spatiotemporal = "off",
  offset = log(sub$swept_area)
)
full_region <- make_survey_prediction_region(sub, crs = 32633, concavity = 5)
full_grid <- sdmTMB::replicate_df(
  make_prediction_grid(full_region, resolution = 20), "year", unique(sub$year)
)
# swept_area is in square nautical miles, so cell area must be too: convert
# the 20 km grid resolution to nautical miles before squaring.
cell_area <- (20 / 1.852)^2
compiled_index <- sdmTMB::get_index_split(
  compiled_model, newdata = full_grid, area = cell_area, nsplit = 3, silent = TRUE
)

# EcoS with the compiled settings, WinterS with its own mesh and fields
sc <- compare_surveys(
  data = sub, compiled_index = compiled_index, compiled_model = compiled_model,
  surveys = list(
    EcoS = "EcoS",
    WinterS = survey_spec("WinterS", mesh_cutoff = 40, spatial = "off")
  ),
  newdata = full_grid, object_crs = 32633,
  mesh_cutoff = 50, offset_col = "swept_area", area = cell_area
)
#> ⠙ Fitting and predicting 2 survey models
#> ✔ Fitting and predicting 2 survey models [7ms]
#> 
sc$settings
#> # A tibble: 2 × 16
#>   survey  series  n_hauls n_years first_year last_year formula            family
#>   <chr>   <chr>     <int>   <int>      <int>     <int> <chr>              <chr> 
#> 1 EcoS    EcoS       1056       6       2020      2025 biomass ~ 0 + as.… tweed…
#> 2 WinterS WinterS    2016       6       2020      2025 biomass ~ 0 + as.… tweed…
#> # ℹ 8 more variables: spatiotemporal <chr>, extra_time <chr>,
#> #   mesh_cutoff <dbl>, mesh_vertices <int>, min_depth <dbl>, max_depth <dbl>,
#> #   grid_cells <int>, shared_cells <int>
plot(sc)

plot(sc, type = "same_area")

plot(sc, type = "scatter")

plot(sc, type = "spatial", fill_prefix = "Biomass density")

# }

Usage

Arguments

Value

See also

Author

Examples