Leave-one-out influence of surveys, cruises or hauls on an sdmTMB index

Refits a fitted sdmTMB::sdmTMB() model once for each identifier in ids with all rows of that identifier removed, and compares the resulting index with the index of a refit to all data. The identifier can be a survey series (which survey carries the signal of a compiled index?), a cruise, or a single haul (which hauls drive an annual index value?).

Each refit reuses the full model's mesh and every setting of its original call (as sdmTMB::update.sdmTMB() does; formula, family, random fields, offset, weights, priors, ...), so only the data change. By default the full model's parameter estimates are used as starting values (falling back to a cold start if the warm-started fit fails).

Removing an identifier can leave years without data, for example when a survey series was the only one sampled in some years. These years are kept in the refit as extra_time, so that the random fields stay defined and the warm start remains valid, and their index values are reported as NA. A formula with a factor year effect cannot be estimated for such years; use a smooth or random-field-only temporal structure, or expect those refits to fail. Likewise, if newdata predicts for a level of a factor in the formula (e.g. cruiseseries = "SeriesA"), that level cannot be left out.

Usage

sdmTMB_leave_out(
  object,
  ids = NULL,
  id_col,
  newdata,
  area = 1,
  index = c("plugin", "split"),
  nsplit = NULL,
  bias_correct = FALSE,
  level = 0.95,
  warm_start = TRUE,
  n_workers = 1,
  future_seed = TRUE,
  verbose = TRUE
)

Arguments

object: A fitted sdmTMB::sdmTMB() model. Its offset, if any, is taken from object$offset.
ids: Values of object$data[[id_col]] to leave out one at a time. If NULL (default), every unique value.
id_col: Name of the identifier column in object$data, e.g. "cruiseseries", a cruise column or a haul identifier.
newdata: Prediction grid with the predictors of the fitted model and its time column.
area: Grid cell area: a single value, a column name in newdata, or a numeric vector with one value per row of newdata.
index: Character. "plugin" (default) sums the response-scale grid predictions times area, which is fast and suits screening many identifiers (e.g. hauls); "split" uses sdmTMB::get_index_split() and returns standard errors and confidence intervals, which suits a few identifiers (e.g. survey series).
nsplit, bias_correct, level: Passed to sdmTMB::get_index_split() when index = "split".
warm_start: Logical. Use the full model's estimates as starting values? Default TRUE.
n_workers: Integer. Number of parallel workers. 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.
verbose: Logical. Show progress messages (via cli when installed)? Default TRUE.

Value

An object of class sdmTMB_leave_out, a list with

full: the index of the refit to all data, by year;
index: the index of each leave-out refit by year, with columns id, year, est (and lwr, upr, se for index = "split"), ratio (est divided by the full index) and removed (TRUE for years with no data left after removing id, where est is NA);
influence: one row per identifier with the number of rows removed, the years it has data in (id_years), haul_share (its share of the rows in those years), influence (one minus the geometric mean ratio over those years, so positive when the identifier raises the index), leverage (influence / haul_share), max_abs_change and year_max_change (the largest absolute relative change over all years) and all_ok (whether the refit passed sdmTMB::sanity()).

Influence and leverage

A large influence on its own mostly flags the main survey of a year: if a cruise supplied 60 % of a year's hauls, removing it is expected to move the index. leverage scales the influence by the identifier's share of the hauls, so a value well above 1 marks an identifier that moves the index far more than its sample size warrants, typically a one-off or auxiliary cruise that sampled a high-density area few other cruises cover. A practical screen is influence >= 0.1 & leverage >= 2; see plot(x, type = "leverage"). Screen candidates cheaply first with density_share().

Author

Mikko Vihtakari

Examples

# \donttest{
data(sebastes)
sub <- subset(sebastes, cruiseseries %in% c("EcoS", "WinterS", "EggaN") &
  year >= 2018)
mesh <- sdmTMB::make_mesh(sub, c("X", "Y"), cutoff = 60)
fit <- sdmTMB::sdmTMB(
  biomass ~ 1, data = sub, time = "year", mesh = mesh,
  family = sdmTMB::tweedie(), spatiotemporal = "ar1",
  offset = log(sub$swept_area)
)
region <- make_survey_prediction_region(sub, crs = 32633)
grid <- sdmTMB::replicate_df(
  make_prediction_grid(region, resolution = 30), "year", unique(sub$year)
)

lo <- sdmTMB_leave_out(
  fit, id_col = "cruiseseries", newdata = grid, area = (30 / 1.852)^2
)
#> ⠙ Refitting the full model
#> ✔ Refitting the full model [4ms]
#> 
#> ⠙ Refitting without each of 3 identifiers
#> ✔ Refitting without each of 3 identifiers [4ms]
#> 
lo$influence
#> # A tibble: 3 × 9
#>   id      n_removed id_years  haul_share influence leverage max_abs_change
#>   <chr>       <int> <list>         <dbl>     <dbl>    <dbl>          <dbl>
#> 1 WinterS      2677 <int [8]>      0.601     0.661    1.10           0.763
#> 2 EcoS         1379 <int [8]>      0.310     0.384    1.24           0.513
#> 3 EggaN         397 <int [4]>      0.160     0.119    0.741          0.324
#> # ℹ 2 more variables: year_max_change <int>, all_ok <lgl>
plot(lo)

plot(lo, type = "ratio")

# }