diff --git a/.Rbuildignore b/.Rbuildignore
index 48e8a6d5..87af7aa4 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -10,4 +10,5 @@
^_pkgdown\.yml$
^docs$
^AGENTS\.md$
+^vignettes/articles$
^\.git$
diff --git a/DESCRIPTION b/DESCRIPTION
index 39ddbba0..a99cc27e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -29,16 +29,16 @@ Imports:
cli,
htmltools,
htmlwidgets,
+ knitr,
+ rmarkdown,
+ flexdashboard,
plotly
Encoding: UTF-8
Roxygen: list(markdown = TRUE)
LazyData: true
Suggests:
pkgload,
- knitr,
- rmarkdown,
jsonlite,
- flexdashboard,
testthat (>= 3.0.0),
withr
Config/testthat/edition: 3
diff --git a/NAMESPACE b/NAMESPACE
index e69a1846..6b6f5716 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -14,6 +14,7 @@ S3method(update,tam_fit)
S3method(update,tam_ref)
S3method(vis_tam,default)
S3method(vis_tam,tam_retro)
+export(add_osa_residuals)
export(add_plot_buttons)
export(ar1)
export(as_tam_data)
@@ -37,6 +38,12 @@ export(fit_retro)
export(fit_tam)
export(for_fleet)
export(for_survey)
+export(get_f)
+export(get_fbar)
+export(get_m)
+export(get_mbar)
+export(get_recruitment)
+export(get_ssb)
export(iid)
export(logistic)
export(make_par)
@@ -47,6 +54,7 @@ export(plot_age_year)
export(plot_bubbles)
export(plot_heatmap)
export(plot_matrix)
+export(plot_model)
export(plot_obs_pred)
export(plot_par)
export(plot_process_deviations)
@@ -84,10 +92,12 @@ importFrom(Matrix,expand)
importFrom(Matrix,solve)
importFrom(cli,cli_abort)
importFrom(cli,cli_warn)
+importFrom(flexdashboard,flex_dashboard)
importFrom(furrr,furrr_options)
importFrom(furrr,future_map)
importFrom(future,multisession)
importFrom(future,plan)
+importFrom(knitr,knit_expand)
importFrom(progressr,progressor)
importFrom(progressr,with_progress)
importFrom(stats,aggregate)
diff --git a/NEWS.md b/NEWS.md
index 16263a2f..6729f5b8 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,32 @@
# tinyAM 0.1.1.9001
+* Add `plot_model()` and a Model structure dashboard page showing biological
+ connections, process switches, formulas, initialization and input predictors.
+ Compare stored specifications with an independent model selector; link
+ components to maintained equations without calculating or changing fits.
+
+* Add Auto / Standardized / OSA controls to public residual plots and dashboards.
+ Auto selects retained OSAs separately for each model and panel, preserving
+ partial gaps and models without OSA calculations. Labels identify the method;
+ changing plot controls never calculates residuals.
+
+* Select any supported OSA algorithm through `method` and retain diagnostics
+ with `add_osa_residuals()`, without refitting. Explicit computation requests
+ refresh catch and survey residuals together; cached extraction stays cheap.
+ Retain method/settings/coverage metadata through tidying and model comparisons.
+
+* Add `get_ssb()`, `get_recruitment()`, `get_fbar()`, `get_mbar()`, `get_f()`
+ and `get_m()` for common population tables, retaining uncertainty, units,
+ projections and comparison labels without calculating diagnostics.
+
+* Print a concise input-validation summary with `check_tam_data()`; internal
+ preparation remains quiet. Install dashboard R dependencies automatically and
+ check for Pandoc before rendering.
+
+* Preserve smooth warm-start functions across equivalent mgcv basis signs and
+ rotations, while rejecting changes to predictor grids or Gaussian penalties.
+ This fixes platform-dependent retrospective warm-start checks.
+
* Read assessment examples from pinned external SAIO snapshots. Separate
curated and preliminary translations, dashboards, caches and comparison outputs.
diff --git a/R/data.R b/R/data.R
index 75960950..11ae73a9 100644
--- a/R/data.R
+++ b/R/data.R
@@ -181,7 +181,7 @@
#' check_tam_data(data)
#' @export
as_tam_data <- function(data) {
- check_tam_data(data)
+ check_tam_data(data, quiet = TRUE)
class(data) <- unique(c("tam_data", class(data)))
data
}
@@ -292,20 +292,26 @@ print.tam_data <- function(x, n = 5, ...) {
#' requires nonempty labels on every row. Units never change input values.
#' See [as_tam_data()] for units of derived quantities.
#'
-#' On failure the function aborts with a `cli` error. On success, it returns
-#' `TRUE` (invisibly).
+#' On failure the function aborts with a `cli` error. On success, it prints a
+#' concise summary and returns `TRUE` invisibly. Use `quiet = TRUE` to suppress
+#' the success summary; validation errors are always reported.
#'
#' @param data A named list expected to contain data.frames `catch`, `index`,
#' `weight`, and `maturity` with the columns described above.
+#' @param quiet Logical. Suppress the success summary? Defaults to `FALSE`.
#'
#' @return `TRUE` (invisibly) if validation passes; otherwise aborts.
#' @examples
#' data(cod_obs)
-#' check_obs(cod_obs) # returns TRUE (invisibly) if valid
+#' check_tam_data(cod_obs)
+#' check_tam_data(cod_obs, quiet = TRUE)
#'
#' @importFrom cli cli_abort cli_warn
#' @export
-check_tam_data <- function(data) {
+check_tam_data <- function(data, quiet = FALSE) {
+ if (!is.logical(quiet) || length(quiet) != 1L || is.na(quiet)) {
+ cli::cli_abort("{.arg quiet} must be TRUE or FALSE.")
+ }
obs <- data
if (!is.list(obs)) {
cli::cli_abort(c(
@@ -478,6 +484,14 @@ check_tam_data <- function(data) {
}
.data_units(obs)
+ if (!quiet) {
+ n_fleets <- if ("fleet" %in% names(obs$catch)) length(unique(obs$catch$fleet)) else 1L
+ n_surveys <- length(unique(obs$index$survey))
+ cli::cli_inform(c(
+ "v" = "Input data passed validation.",
+ "i" = "Years {years[1]}-{years[length(years)]}; ages {ages[1]}-{ages[length(ages)]}; fleets: {n_fleets}; surveys: {n_surveys}."
+ ))
+ }
invisible(TRUE)
}
@@ -513,10 +527,11 @@ check_tam_data <- function(data) {
#'
#' A convenience wrapper for [check_tam_data()].
#' @param obs A named list of tinyAM input tables.
+#' @inheritParams check_tam_data
#' @inherit check_tam_data details
#' @return `TRUE` (invisibly) if validation passes; otherwise aborts.
#' @examples
#' check_tam_data(cod_obs)
#' @export
-check_obs <- function(obs) check_tam_data(obs)
+check_obs <- function(obs, quiet = FALSE) check_tam_data(obs, quiet = quiet)
diff --git a/R/fit.R b/R/fit.R
index 2e5d8a6f..13c43dc9 100644
--- a/R/fit.R
+++ b/R/fit.R
@@ -105,11 +105,14 @@
start$log_f <- start$log_f[, groups[keep], drop = FALSE]
}
}
+ smooth_start <- .smooth_warm_starts(start, source$dat, dat)
+ start <- smooth_start$parameters
for (component in c("N", "F", "M", "rec", "q")) {
name <- switch(component, N = "mu_n", F = "log_mu_f", M = "mu_m", rec = "rec_beta",
q=if (dat$index_settings$q_link == "logit") "logit_q" else "log_q")
if (is.null(start[[name]])) next
old <- if (component == "rec") source$dat$rec$matrix else source$dat[[paste0(component, "_modmat")]]
+ if (component %in% c("F", "q")) old <- smooth_start$designs[[component]]
new <- if (component == "rec") dat$rec$matrix else dat[[paste0(component, "_modmat")]]
if (!is.matrix(new)) next
if (component %in% c("F","q")) {
@@ -138,15 +141,6 @@
cli::cli_inform("{label} mean starts were omitted because the formula basis changed.")
}
}
- for (component in c("q","F")) {
- field <- if (component == "q") "q_terms" else "F_terms"
- for (term in Filter(function(x)x$type == "smooth",dat[[field]])) {
- previous <- source$dat[[field]][[term$id]]
- if (is.null(previous) || !.smooth_basis_compatible(previous,term)) {
- start[c(term$parameter,term$sd_parameter)] <- NULL
- }
- }
- }
}
if (!is.list(start) || (length(start) && (is.null(names(start)) ||
anyNA(names(start)) || any(!nzchar(names(start))) || anyDuplicated(names(start))))) {
@@ -222,9 +216,9 @@
#' @param interval Level in `(0, 1)` to use to generate confidence intervals,
#' where applicable; default `0.95.`
#' @param silent Logical; if `TRUE`, disables tracing information.
-#' @param add_osa_res Logical; add one-step-ahead residuals? Hard wired to
-#' apply the `"oneStepGaussianOffMode"` method. See [RTMB::oneStepPredict()]
-#' for details.
+#' @param add_osa_res Logical; add one-step-ahead residuals using
+#' `"oneStepGaussianOffMode"`? To select another method after fitting and
+#' retain its results, use [add_osa_residuals()].
#' @param start_par Optional fitted `tam_fit` or list of starting parameter values, matching the
#' structure returned by [make_par()]. Useful for warm starts or retrospective
#' runs where parameter estimates from a previous fit are reused to improve
diff --git a/R/getters.R b/R/getters.R
new file mode 100644
index 00000000..5e32a3a9
--- /dev/null
+++ b/R/getters.R
@@ -0,0 +1,104 @@
+#' Extract common population estimates
+#'
+#' Return a ready-to-plot table of spawning biomass, recruitment or mortality,
+#' without navigating the nested tidy output. These functions select the same
+#' estimates, uncertainty and units used in [tidy_tam()] and the dashboard.
+#'
+#' @param ... One or more fitted `tam_fit` models or `tam_ref` reporting
+#' references, or one named list of these objects.
+#' @param model_list A named list of models, instead of `...`.
+#' @param interval Confidence level in `(0, 1)`, default `0.95`. Fitted summaries
+#' are rebuilt when their cached interval differs. Reference uncertainty is
+#' retained as supplied, because it cannot be recalculated from a reporting
+#' reference.
+#' @details
+#' `get_ssb()` and `get_recruitment()` return annual estimates. `get_fbar()`
+#' and `get_mbar()` return the model's abundance-weighted means over the ages
+#' chosen by `F_settings$mean_ages` and `M_settings$mean_ages`, respectively.
+#' `get_f()` and `get_m()` return year-by-age surfaces; F is total fishing
+#' mortality, summed over fleets. They do not return process deviations.
+#'
+#' A single fit returns its table without a model column. Multiple fits and
+#' named-list inputs add `model`, using the supplied names. Missing reference
+#' quantities are omitted rather than reconstructed. Explicit reference
+#' `comparison_scales` follow [tidy_tam()]'s conversion rules; otherwise known
+#' incompatible units must be converted before combining models.
+#'
+#' These are lightweight selectors. They never refit models or calculate
+#' observation OSA residuals or joint-sample process diagnostics. Supplied
+#' objects are not modified.
+#'
+#' @return A data frame retaining available keys, `est`, confidence limits,
+#' `se`, `se_scale`, `unit` and `is_proj`. Available metadata, such as SSB
+#' timing, is retained. Annual summaries have `year`; surfaces also have `age`.
+#' @examples
+#' \dontrun{
+#' get_ssb(fit)
+#' get_recruitment(fit, interval = 0.8)
+#' get_fbar(Base = fit, Alternative = other_fit)
+#' get_f(model_list = list(Base = fit, Alternative = other_fit))
+#' }
+#' @seealso [tidy_pop()], [tidy_tam()], [plot_trend()], [plot_matrix()]
+#' @name get_population
+NULL
+
+#' @rdname get_population
+#' @export
+get_ssb <- function(..., model_list = NULL, interval = 0.95) {
+ .get_population_quantity("ssb", list(...), as.list(substitute(list(...)))[-1L],
+ model_list, interval)
+}
+
+#' @rdname get_population
+#' @export
+get_recruitment <- function(..., model_list = NULL, interval = 0.95) {
+ .get_population_quantity("recruitment", list(...), as.list(substitute(list(...)))[-1L],
+ model_list, interval)
+}
+
+#' @rdname get_population
+#' @export
+get_fbar <- function(..., model_list = NULL, interval = 0.95) {
+ .get_population_quantity("F_bar", list(...), as.list(substitute(list(...)))[-1L],
+ model_list, interval)
+}
+
+#' @rdname get_population
+#' @export
+get_mbar <- function(..., model_list = NULL, interval = 0.95) {
+ .get_population_quantity("M_bar", list(...), as.list(substitute(list(...)))[-1L],
+ model_list, interval)
+}
+
+#' @rdname get_population
+#' @export
+get_f <- function(..., model_list = NULL, interval = 0.95) {
+ .get_population_quantity("F", list(...), as.list(substitute(list(...)))[-1L],
+ model_list, interval)
+}
+
+#' @rdname get_population
+#' @export
+get_m <- function(..., model_list = NULL, interval = 0.95) {
+ .get_population_quantity("M", list(...), as.list(substitute(list(...)))[-1L],
+ model_list, interval)
+}
+
+.get_population_quantity <- function(quantity, dots, expressions, model_list, interval) {
+ if (!is.numeric(interval) || length(interval) != 1L || !is.finite(interval) ||
+ interval <= 0 || interval >= 1) {
+ cli::cli_abort("{.arg interval} must be one number strictly between zero and one.")
+ }
+ info <- .dots_or_list(dots, expressions, model_list, "model_list")
+ named_list <- length(dots) == 1L && is.list(dots[[1L]]) &&
+ !.is_reporting_fit(dots[[1L]])
+ add_label <- !info$using_dots || named_list || length(info$fits) > 1L
+ populations <- .population_tables(info$fits, interval)
+ selected <- lapply(populations, function(x) x[intersect(names(x), quantity)])
+ tab <- stack_nested(selected, label = if (add_label) "model" else NULL)[[quantity]]
+ if (!is.data.frame(tab)) {
+ cli::cli_abort("The requested quantity {.val {quantity}} is unavailable in the supplied models.")
+ }
+ .unit_label(tab, "")
+ tab
+}
diff --git a/R/model-diagram.R b/R/model-diagram.R
new file mode 100644
index 00000000..039e876b
--- /dev/null
+++ b/R/model-diagram.R
@@ -0,0 +1,282 @@
+#' View the structure of fitted assessment models
+#'
+#' See how recruitment, survival, mortality and observations are connected,
+#' and which components have annual process variation. The diagram and its
+#' accompanying cards use the completed settings stored in each fitted model.
+#'
+#' @param ... One or more `tam_fit` models or `tam_ref` reporting references,
+#' or one named list of these objects.
+#' @param model_list A named list of models, instead of `...`.
+#' @details
+#' Multiple models have a selector, initially showing the first tinyAM fit.
+#' A reporting reference displays an unavailable message because it does not
+#' contain a fitted tinyAM specification. The same viewer appears on the
+#' dashboard's Model structure page.
+#'
+#' An annual process set to `"off"` is labeled deterministic: it removes
+#' process deviations, not the underlying recruitment, survival or mortality.
+#' Initial states are shown separately from annual processes. Stock-recruit
+#' relationships add a parent-SSB feedback arrow. Formula cards retain the
+#' actual formulas, fleets, surveys, state blocks and projection controls.
+#' Short formula labels appear in the diagram; dashed input arrows identify
+#' the catch, index or population table supplying active predictors.
+#'
+#' This is a model overview, not a diagnostic of identifiability or numerical
+#' convergence. It reads stored settings without evaluating the fitted
+#' objective, sampling uncertainty, calculating residuals or changing the fit.
+#' No additional diagramming package is required.
+#'
+#' @return A browsable HTML container with a responsive SVG diagram and model
+#' settings. Display it directly, or save it with [htmltools::save_html()].
+#' @examples
+#' \dontrun{
+#' plot_model(fit)
+#' plot_model(Base = fit, Alternative = other_fit)
+#' }
+#' @seealso [vis_tam()], [fit_tam()], [tinyAM-model]
+#' @export
+plot_model <- function(..., model_list = NULL) {
+ info <- .dots_or_list(list(...), as.list(substitute(list(...)))[-1L],
+ model_list, "model_list")
+ fits <- info$fits
+ first <- which(vapply(fits, inherits, logical(1), "tam_fit"))
+ first <- if (length(first)) first[1L] else 1L
+ panes <- lapply(seq_along(fits), function(i) {
+ content <- if (inherits(fits[[i]], "tam_fit")) {
+ .model_structure_view(fits[[i]]$dat, names(fits)[i])
+ } else {
+ htmltools::tags$p(role = "status", class = "tam-model-unavailable",
+ "Model structure is unavailable for a reporting reference (tam_ref).")
+ }
+ htmltools::tags$div(`data-model-pane` = as.character(i),
+ style = if (i != first) "display:none" else NULL, content)
+ })
+ selector <- if (length(fits) > 1L) htmltools::tags$label(
+ "Model ", htmltools::tags$select(`aria-label` = "Model structure model",
+ onchange = paste0("var c=this.closest('.tam-model-structure');",
+ "var v=this.value;c.querySelectorAll('[data-model-pane]').forEach(",
+ "function(p){p.style.display=p.getAttribute('data-model-pane')===v?'':'none';});"),
+ lapply(seq_along(fits), function(i) htmltools::tags$option(
+ value = as.character(i), selected = if (i == first) "selected" else NULL,
+ names(fits)[i]))))
+ htmltools::browsable(htmltools::tags$div(class = "tam-model-structure",
+ htmltools::tags$style(htmltools::HTML(paste0(
+ ".tam-model-structure{font-family:Arial,sans-serif;color:#1f2933;background:#fff;padding:8px;overflow:auto}",
+ ".tam-model-structure select{margin:0 0 12px 8px;padding:5px;max-width:80%}",
+ ".tam-model-structure svg{width:100%;height:auto;min-width:630px;display:block}",
+ ".tam-model-cards{display:grid;grid-template-columns:repeat(auto-fit,minmax(235px,1fr));gap:12px}",
+ ".tam-model-card{border:1px solid #d6dce0;border-top:3px solid #007187;padding:10px}",
+ ".tam-model-card h4{font-size:16px;color:#003e56;margin:0 0 6px}",
+ ".tam-model-card p{font-size:13px;margin:4px 0;overflow-wrap:anywhere}",
+ ".tam-model-card code{font-size:12px;white-space:pre-wrap}",
+ ".tam-model-legend,.tam-model-note{font-size:13px;margin:8px 0 12px}",
+ ".tam-model-unavailable{padding:25px;border:1px solid #d6dce0}"))),
+ selector, panes))
+}
+
+.model_structure_spec <- function(dat) {
+ form <- function(x) if (is.null(x)) "None" else paste(deparse(x, width.cutoff = 90L), collapse = " ")
+ process <- function(x) if (identical(x, "off")) "Deterministic (process off)" else
+ paste("Stochastic process:", if (is.null(x)) "unspecified" else toupper(x))
+ block_text <- function(x) {
+ if (is.null(x)) return("Not specified")
+ if (is.list(x)) return(paste(vapply(names(x), function(nm)
+ paste0(nm, ": ", block_text(x[[nm]])), character(1)), collapse = "; "))
+ if (is.factor(x)) return(paste(levels(x), collapse = ", "))
+ paste(x, collapse = ", ")
+ }
+ fleets <- if (isTRUE(dat$has_fleet)) dat$fleets else "Aggregate fishery"
+ surveys <- unique(as.character(dat$obs$index$survey))
+ curve <- dat$rec$curve
+ predictors <- function(forms, table, data, exclude = character()) {
+ forms <- Filter(function(x) inherits(x, "formula"), forms)
+ variables <- unique(unlist(lapply(forms, .formula_input_variables), use.names = FALSE))
+ # Population year/age coordinates exist even without pop_covariates.
+ if (table == "pop_covariates") exclude <- c(exclude, "year", "age")
+ variables <- intersect(setdiff(variables, c("obs", exclude)), names(data))
+ if (length(variables)) paste0(table, ": ", paste(variables, collapse = ", ")) else ""
+ }
+ correlation <- if (inherits(dat$F_settings$cor_form, "formula"))
+ list(dat$F_settings$cor_form) else dat$F_settings$cor_form
+ inputs <- c(
+ recruitment = predictors(list(dat$rec_settings$mu_form), "pop_covariates",
+ dat$pop_data, exclude = "ssb"),
+ population = predictors(list(dat$N_settings$mu_form,
+ if (dat$N_settings$process != "off") dat$N_settings$sd_form),
+ "pop_covariates", dat$pop_data),
+ fishing = predictors(c(list(dat$F_settings$mu_form),
+ if (dat$F_settings$process != "off") list(dat$F_settings$sd_form),
+ if (dat$F_settings$process %in% c("ar1", "cor_rw")) correlation),
+ "catch", dat$F_data),
+ natural = predictors(list(dat$M_settings$mu_form, dat$M_settings$mu_supplied,
+ if (dat$M_settings$process != "off") dat$M_settings$sd_form),
+ "pop_covariates", dat$pop_data),
+ catch = predictors(list(dat$catch_settings$sd_form, dat$catch_settings$sd_supplied),
+ "catch", dat$obs$catch),
+ survey = predictors(list(dat$index_settings$q_form, dat$index_settings$sd_form,
+ dat$index_settings$sd_supplied), "index", dat$obs$index, exclude = "N"))
+ formula_labels <- c(initial = paste0("N0: ", dat$N_settings$init,
+ "; R0: ", dat$rec_settings$init),
+ recruitment = if (is.null(dat$rec_settings$mu_form)) "Mean: previous year" else
+ paste("Mean", form(dat$rec_settings$mu_form)),
+ population = if (is.null(dat$N_settings$mu_form)) "Mean: survivors" else
+ paste("Mean", form(dat$N_settings$mu_form)),
+ fishing = if (is.null(dat$F_settings$mu_form)) "Mean: previous year" else
+ paste("Mean", form(dat$F_settings$mu_form)),
+ natural = if (is.null(dat$M_settings$mu_form)) paste("Baseline", form(dat$M_settings$mu_supplied)) else
+ paste("Mean effects", form(dat$M_settings$mu_form)),
+ catch = paste("SD", form(dat$catch_settings$sd_form)),
+ survey = paste("q", form(dat$index_settings$q_form)),
+ ssb = "SSB = sum of mature biomass")
+ card <- function(title, lines, help) list(title = title, lines = lines, help = help)
+ cards <- list(
+ initial = card("Initial states", c(
+ paste("Older-age initialization:", dat$N_settings$init),
+ paste("Recruitment boundary:", dat$rec_settings$init),
+ paste("Initial year:", min(dat$years)),
+ "Initial parameters are separate from annual process deviations."), "recruitment_formulas"),
+ recruitment = card("Recruitment", c(process(dat$rec_settings$process),
+ paste("Mean:", form(dat$rec_settings$mu_form)),
+ if (!is.null(curve)) paste("Parent-SSB relationship:", curve$type, "; lag", curve$lag),
+ "Covariates: pop_covariates."), "recruitment_formulas"),
+ population = card("Abundance and cohort survival", c(process(dat$N_settings$process),
+ paste("Mean:", form(dat$N_settings$mu_form)),
+ if (!identical(dat$N_settings$process, "off")) paste("SD:", form(dat$N_settings$sd_form)),
+ paste("Modeled ages:", paste(range(dat$ages), collapse = " to "), "; oldest age is the plus group."),
+ "Apparent gains/losses do not identify migration separately from mortality."), "tinyAM-model"),
+ fishing = card("Fishing mortality", c(process(dat$F_settings$process),
+ paste("Mean:", form(dat$F_settings$mu_form)),
+ if (!identical(dat$F_settings$process, "off")) paste("SD:", form(dat$F_settings$sd_form)),
+ if (dat$F_settings$process %in% c("ar1", "cor_rw")) paste("Correlation groups:", form(dat$F_settings$cor_form)),
+ paste("Fleets:", paste(fleets, collapse = ", ")),
+ paste("Age blocks:", block_text(if (isTRUE(dat$has_fleet)) dat$F_blocks else dat$F_settings$age_blocks)),
+ if (any(dat$F_settings$state_index == 0L)) "F is exactly zero outside the declared fishing-age blocks.",
+ "Total F is the sum of fleet F. Shared age blocks share actual F."), "for_fleet"),
+ natural = card("Natural mortality", c(process(dat$M_settings$process),
+ paste("Supplied mean:", form(dat$M_settings$mu_supplied)),
+ paste("Mean effects:", form(dat$M_settings$mu_form)),
+ if (!identical(dat$M_settings$process, "off")) paste("SD:", form(dat$M_settings$sd_form)),
+ paste("Age blocks:", block_text(dat$M_settings$age_blocks)),
+ if (!identical(dat$M_settings$process, "off")) paste("First process year:", dat$M_settings$first_dev_year),
+ "Shared age blocks share actual M; other cells retain their mean."), "tinyAM-model"),
+ catch = card("Catch observations", c("Lognormal numbers at age, conditional on population and F.",
+ paste("Observation SD:", form(dat$catch_settings$sd_form)),
+ paste("Supplied SD:", form(dat$catch_settings$sd_supplied)),
+ "Zero and missing observations are excluded from the lognormal likelihood."), "fit_tam"),
+ survey = card("Survey indices and catchability", c(
+ paste("Surveys:", paste(surveys, collapse = ", ")),
+ paste("q formula:", form(dat$index_settings$q_form)),
+ paste("q link:", dat$index_settings$q_link),
+ paste("Gaussian q terms:", {
+ types <- vapply(dat$q_terms, `[[`, character(1), "type")
+ types <- unique(types[types %in% c("iid", "rw", "ar1", "smooth")])
+ if (length(types)) paste(types, collapse = ", ") else "None"
+ }),
+ paste("Observation SD:", form(dat$index_settings$sd_form)),
+ paste("Supplied SD:", form(dat$index_settings$sd_supplied)),
+ "Predictions use each survey's sampling time; q may contain Gaussian effects."), "for_fleet"),
+ ssb = card("Spawning biomass and projections", c(
+ "SSB uses stock weights, maturity and spawning timing, retaining plus-group biology.",
+ if (is.null(dat$proj_settings) || !isTRUE(dat$proj_settings$n_proj > 0)) "Projections: none." else
+ paste("Projection settings:", form(dat$proj_settings))), "tinyAM-model"))
+ list(cards = cards, feedback = !is.null(curve),
+ feedback_label = if (!is.null(curve)) paste0("Parent SSB: ", curve$type,
+ ", lag ", curve$lag) else "", inputs = inputs, formula_labels = formula_labels,
+ statuses = c(recruitment = dat$rec_settings$process,
+ population = dat$N_settings$process, fishing = dat$F_settings$process,
+ natural = dat$M_settings$process),
+ fleets = fleets, surveys = surveys)
+}
+
+.model_structure_view <- function(dat, label) {
+ spec <- .model_structure_spec(dat)
+ escape <- function(x) as.character(htmltools::htmlEscape(x))
+ shorten <- function(x, width = 34L) if (nchar(x) > width)
+ paste0(substr(x, 1L, width - 3L), "...") else x
+ positions <- list(initial = c(30, 35), recruitment = c(30, 200),
+ fishing = c(30, 370), natural = c(310, 370), population = c(590, 200),
+ catch = c(920, 370), survey = c(920, 35), ssb = c(590, 425))
+ subtitles <- c(initial = "Boundary parameters", recruitment = "New youngest-age fish",
+ fishing = "Total F = sum of fleet F", natural = "Shared population M",
+ population = "Survivors + recruitment", catch = "Fleet catch predictions",
+ survey = "q, sampling time and abundance", ssb = "Mature stock biomass")
+ titles <- c(initial = "Initial states", recruitment = "Recruitment",
+ fishing = "Fishing mortality", natural = "Natural mortality",
+ population = "Population abundance", catch = "Catch observations",
+ survey = "Survey observations", ssb = "Spawning biomass")
+ arrow <- function(points) {
+ n <- nrow(points)
+ direction <- points[n, ] - points[n - 1L, ]
+ direction <- direction / sqrt(sum(direction^2))
+ tip <- points[n, ]; base <- tip - 9 * direction
+ normal <- c(-direction[2L], direction[1L]) * 4
+ paste0('',
+ '')
+ }
+ edges <- list(rbind(c(270, 90), c(550, 90), c(550, 215), c(590, 215)),
+ rbind(c(270, 255), c(590, 255)),
+ rbind(c(270, 425), c(295, 425), c(295, 335), c(640, 335), c(640, 310)),
+ rbind(c(550, 425), c(570, 425), c(570, 290), c(590, 290)),
+ rbind(c(830, 230), c(875, 230), c(875, 90), c(920, 90)),
+ rbind(c(830, 280), c(875, 280), c(875, 425), c(920, 425)),
+ rbind(c(150, 480), c(150, 548), c(1040, 548), c(1040, 480)),
+ rbind(c(710, 310), c(710, 425)))
+ if (spec$feedback) edges[[length(edges) + 1L]] <-
+ rbind(c(590, 480), c(285, 480), c(285, 275), c(270, 275))
+ nodes <- vapply(names(positions), function(id) {
+ p <- positions[[id]]
+ status <- spec$statuses[id]
+ stochastic <- length(status) && !is.na(status) && status != "off"
+ fill <- if (stochastic) "#d9efed" else "#f0f3f5"
+ annotation <- if (is.na(status)) "" else if (stochastic)
+ paste("Stochastic", toupper(status)) else "Deterministic: process off"
+ title <- paste(spec$cards[[id]]$lines, collapse = "\n")
+ paste0('', escape(title), '',
+ '',
+ '', escape(titles[id]), '',
+ '', escape(subtitles[id]), '',
+ '', escape(annotation), '',
+ '', escape(shorten(spec$formula_labels[id])), '')
+ }, character(1))
+ input_arrows <- vapply(names(spec$inputs)[nzchar(spec$inputs)], function(id) {
+ p <- positions[[id]]
+ paste0('', escape(spec$inputs[id]), '',
+ '',
+ '',
+ escape(shorten(spec$inputs[id], 36L)), '',
+ '',
+ '')
+ }, character(1))
+ svg <- paste0('')
+ cards <- lapply(spec$cards, function(card) htmltools::tags$section(class = "tam-model-card",
+ htmltools::tags$h4(card$title), lapply(card$lines, htmltools::tags$p),
+ htmltools::tags$a(href = paste0("https://paulregular.github.io/tinyAM/reference/", card$help, ".html"),
+ "Equations and settings")))
+ htmltools::tagList(htmltools::tags$p(class = "tam-model-note",
+ "Stored specification for ", htmltools::tags$strong(label),
+ ". Arrows show biological and observation relationships, not a computation order."),
+ htmltools::HTML(svg), htmltools::tags$p(class = "tam-model-legend",
+ "Teal: annual stochastic process. Gray: no annual process deviations. ",
+ "Dashed arrows: active input-table predictors. Process off retains the biological component. ",
+ "Initial parameters and observations are separate."),
+ htmltools::tags$details(htmltools::tags$summary("Formulas and settings"),
+ htmltools::tags$div(class = "tam-model-cards", cards)))
+}
diff --git a/R/osa.R b/R/osa.R
new file mode 100644
index 00000000..481ceed4
--- /dev/null
+++ b/R/osa.R
@@ -0,0 +1,197 @@
+#' Calculate and retain one-step-ahead residuals
+#'
+#' Calculate observation diagnostics at fitted parameter values and return a fit
+#' containing the results. Use the returned object in [vis_tam()] or [residuals()].
+#' No model is refitted, and the supplied fit is not modified.
+#'
+#' @param fit A fitted `tam_fit`.
+#' @param method One algorithm supported by [RTMB::oneStepPredict()]:
+#' `"oneStepGaussianOffMode"` (default), `"oneStepGaussian"`,
+#' `"oneStepGeneric"`, `"fullGaussian"`, or `"cdf"`.
+#' @param ... Named computation controls for [RTMB::oneStepPredict()], such as
+#' `trace`, `parallel`, `subset`, `conditional`, or `splineApprox`.
+#' `subset` and `conditional` index included log observations in likelihood
+#' order. Their unrequested rows retain missing residuals.
+#' @details The algorithms use different predictive approximations and can
+#' differ numerically. `oneStepGeneric` uses numerical integration and can be
+#' expensive. The existing off-mode Gaussian calculation remains the default.
+#' Log observations have a continuous Gaussian likelihood, so `discrete = FALSE`
+#' is supplied unless explicitly overridden. Method-dependent upstream defaults,
+#' including `reverse`, are retained. When SAM is also loaded, parallel OSA
+#' calculations use the existing sequential safeguard.
+#'
+#' Temporary RTMB objectives protect the fitted parameters, uncertainty and
+#' objective state. Only OSA columns and metadata are added to a copy of the
+#' existing observation tables; row order, confidence intervals, units and custom
+#' columns are preserved. Metadata in `fit$osa` records the method, requested and
+#' effective settings, and observation coverage. Failed requested calculations
+#' remain missing and are recorded separately from unrequested observations.
+#' A calculation error returns no partially updated fit.
+#' @return A `tam_fit` with `osa_res` and `osa_method` columns in `$obs_pred`,
+#' and calculation metadata in `$osa` and the observation-table attributes.
+#' @examples
+#' \dontrun{
+#' fit_osa <- add_osa_residuals(fit, method = "oneStepGeneric", trace = FALSE)
+#' residuals(fit_osa, type = "osa")
+#' vis_tam(fit_osa)
+#' }
+#' @seealso [residuals.tam_fit()], [tidy_obs_pred()]
+#' @export
+add_osa_residuals <- function(fit, method = "oneStepGaussianOffMode", ...) {
+ fit <- .require_tam_fit(fit, arg = "fit")
+ result <- .compute_osa_residuals(fit, method, list(...))
+ tables <- .observation_tables(fit)
+ tables <- .attach_osa_tables(tables, result$tables, result$metadata)
+ out <- fit
+ out$obs_pred <- tables
+ out$osa <- result$metadata
+ out
+}
+
+.osa_method <- function(method, controls = list()) {
+ if ("methods" %in% names(controls)) {
+ cli::cli_abort("Use {.arg method}, not {.arg methods}, to select one OSA algorithm.")
+ }
+ supported <- c("oneStepGaussianOffMode", "oneStepGaussian", "oneStepGeneric",
+ "fullGaussian", "cdf")
+ if (!is.character(method) || length(method) != 1L || is.na(method) ||
+ !method %in% supported) {
+ cli::cli_abort("{.arg method} must be one of: {paste(supported, collapse = ', ')}.")
+ }
+ if (length(controls) && (is.null(names(controls)) || anyNA(names(controls)) ||
+ any(!nzchar(names(controls))) || anyDuplicated(names(controls)))) {
+ cli::cli_abort("OSA computation controls must have unique, nonempty names.")
+ }
+ reserved <- intersect(names(controls), c("obj", "method", "observation.name", "data.term.indicator"))
+ if (length(reserved)) {
+ cli::cli_abort("These OSA arguments are supplied by tinyAM: {paste(reserved, collapse = ', ')}.")
+ }
+ method
+}
+
+.osa_indices <- function(x, n, name) {
+ if (is.null(x)) return(integer())
+ if (!is.numeric(x) || anyNA(x) || any(!is.finite(x)) ||
+ any(x < 1L | x > n) || any(x != floor(x)) || anyDuplicated(x)) {
+ cli::cli_abort("OSA {.arg {name}} must contain distinct observation indices between 1 and {n}.")
+ }
+ as.integer(x)
+}
+
+.compute_osa_residuals <- function(fit, method, controls = list()) {
+ method <- .osa_method(method, controls)
+ if (is.null(fit$obj) || is.null(fit$parameter_values)) {
+ cli::cli_abort("OSA calculation requires the fitted objective and stored parameter values.")
+ }
+ if (is.null(fit$dat$obs_map) || is.null(fit$dat$observed)) {
+ cli::cli_abort("OSA calculation requires the stored observations and likelihood mapping.")
+ }
+ requested_settings <- controls
+ n <- length(fit$dat$observed)
+ conditional <- .osa_indices(controls$conditional, n, "conditional")
+ selected <- if (is.null(controls$subset)) setdiff(seq_len(n), conditional) else
+ .osa_indices(controls$subset, n, "subset")
+ if (length(intersect(selected, conditional))) {
+ cli::cli_abort("OSA {.arg subset} and {.arg conditional} must not overlap.")
+ }
+ if (is.null(controls$discrete)) controls$discrete <- FALSE
+ if (isTRUE(controls$parallel) && "stockassessment" %in% loadedNamespaces()) {
+ controls$parallel <- FALSE
+ cli::cli_inform("Calculating one-step residuals sequentially while SAM and RTMB are loaded.")
+ }
+ values <- rep(NA_real_, n)
+ if (length(selected)) {
+ result <- tryCatch({
+ obj <- .fitted_objective(fit)
+ do.call(RTMB::oneStepPredict, c(list(obj = obj, method = method), controls))
+ }, error = function(e) cli::cli_abort(
+ "Could not calculate OSA residuals with {method}: {conditionMessage(e)}", parent = e))
+ if (!is.numeric(result$residual) || length(result$residual) != length(selected)) {
+ cli::cli_abort("OSA results do not match the requested observation subset.")
+ }
+ values[selected] <- result$residual
+ }
+ failed <- selected[!is.finite(values[selected])]
+ if (length(failed)) {
+ values[failed] <- NA_real_
+ cli::cli_warn("OSA method {method} returned unavailable residuals for {length(failed)} requested observations.")
+ }
+ map <- fit$dat$obs_map
+ observed <- map$is_observed
+ coverage <- map[observed, intersect(c("type", "year", "age", "survey", "fleet"), names(map)), drop = FALSE]
+ coverage$observation <- seq_len(n)
+ coverage$requested <- seq_len(n) %in% selected
+ coverage$conditional <- seq_len(n) %in% conditional
+ coverage$available <- is.finite(values)
+ coverage$status <- ifelse(coverage$conditional, "conditional",
+ ifelse(!coverage$requested, "not_requested", ifelse(coverage$available, "available", "failed")))
+ rownames(coverage) <- NULL
+ metadata <- list(method = method, requested_settings = requested_settings,
+ settings = c(list(method = method), controls), coverage = coverage,
+ n_observations = n, n_requested = length(selected), n_available = sum(is.finite(values)),
+ n_failed = length(failed), n_excluded = sum(!observed))
+ all_values <- rep(NA_real_, nrow(map))
+ all_values[observed] <- values
+ tables <- fit$dat$obs[c("catch", "index")]
+ for (name in names(tables)) {
+ tables[[name]]$osa_res <- all_values[map$type == name]
+ tables[[name]]$osa_method <- rep(method, nrow(tables[[name]]))
+ attr(tables[[name]], "osa") <- metadata
+ }
+ attr(tables, "osa") <- metadata
+ list(tables = tables, metadata = metadata)
+}
+
+.observation_tables <- function(fit) {
+ tables <- fit$obs_pred
+ absent <- vapply(c("catch", "index"), function(name) !is.data.frame(tables[[name]]), logical(1))
+ if (any(absent)) {
+ interval <- attr(tables, "interval", exact = TRUE)
+ if (is.null(interval)) interval <- attr(fit$pop, "interval", exact = TRUE)
+ if (is.null(interval)) interval <- .95
+ fallback <- tidy_obs_pred(fit, interval = interval)
+ for (name in names(absent)[absent]) tables[[name]] <- fallback[[name]]
+ if (is.null(attr(tables, "interval", exact = TRUE))) attr(tables, "interval") <- interval
+ }
+ tables
+}
+
+.observation_match <- function(target, source) {
+ keys <- intersect(c("year", "age", "survey", "fleet"), names(source))
+ if (!all(keys %in% names(target)) || !length(keys)) {
+ cli::cli_abort("OSA observation rows must retain the fitted year, age, survey and fleet keys.")
+ }
+ key <- function(d) do.call(paste, c(d[keys], sep = "\r"))
+ source_key <- key(source)
+ rows <- match(key(target), source_key)
+ if (anyDuplicated(source_key) || anyNA(rows)) {
+ cli::cli_abort("OSA observation rows do not match the fitted inputs uniquely.")
+ }
+ rows
+}
+
+.attach_osa_tables <- function(tables, source, metadata = attr(source, "osa", exact = TRUE)) {
+ for (name in intersect(c("catch", "index"), names(tables))) {
+ if (!is.data.frame(tables[[name]]) || !is.data.frame(source[[name]]) ||
+ !"osa_res" %in% names(source[[name]])) next
+ rows <- .observation_match(tables[[name]], source[[name]])
+ tables[[name]]$osa_res <- source[[name]]$osa_res[rows]
+ tables[[name]]$osa_method <- if ("osa_method" %in% names(source[[name]]))
+ as.character(source[[name]]$osa_method[rows]) else rep(
+ if (is.null(metadata$method)) NA_character_ else metadata$method, nrow(tables[[name]]))
+ if (!is.null(metadata)) attr(tables[[name]], "osa") <- metadata
+ }
+ if (!is.null(metadata)) attr(tables, "osa") <- metadata
+ tables
+}
+
+.preserve_observation_extras <- function(tables, cached, metadata = attr(cached, "osa", exact = TRUE)) {
+ for (name in intersect(names(tables), names(cached))) {
+ extras <- setdiff(names(cached[[name]]), names(tables[[name]]))
+ if (length(extras)) {
+ rows <- .observation_match(tables[[name]], cached[[name]])
+ tables[[name]][extras] <- cached[[name]][rows, extras, drop = FALSE]
+ }
+ }
+ .attach_osa_tables(tables, cached, metadata)
+}
diff --git a/R/plot.R b/R/plot.R
index 4835dcc3..fd08da2d 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -79,6 +79,12 @@
#' @param res Optional one-column formula for signed bubble values, e.g.
#' `~est` for process deviations. Omitted values use OSA or standardized
#' observation residuals. Use `y = ~age_block` for shared M state blocks.
+#' @param residual_type Initial observation residual view: `"auto"` (default),
+#' `"standardized"` or `"osa"`. Auto uses available cached OSA residuals
+#' separately for each model in the supplied panel, otherwise standardized
+#' residuals. Partial OSA gaps remain missing. A dropdown switches views;
+#' unavailable OSA values are identified without computing them. A supplied
+#' `res` formula in `plot_bubbles()` bypasses these observation controls.
#' @param ylab Label for the y-axis.
#' @param xlab Label for the x-axis.
#' @param zlab Label for the colorbar in heatmaps.
@@ -411,7 +417,12 @@ plot_obs_pred <- function(
#' @rdname plot_tam
#' @export
-plot_bubbles <- function(data, res = NULL, y = ~age, ylab = "Age", title = NULL, ...) {
+plot_bubbles <- function(data, res = NULL, y = ~age, ylab = "Age", title = NULL,
+ residual_type = c("auto", "standardized", "osa"), ...) {
+ if (is.null(res)) {
+ return(.plot_observation_residuals(data, match.arg(residual_type), "bubbles",
+ list(...), y = y, ylab = ylab, title = title))
+ }
if (!is.null(res)) {
if (!inherits(res, "formula") || length(res) != 2L) cli::cli_abort("res must be a one-sided formula selecting a numeric column.")
values <- stats::model.frame(res, data, na.action = stats::na.pass)
@@ -419,12 +430,6 @@ plot_bubbles <- function(data, res = NULL, y = ~age, ylab = "Age", title = NULL,
data$res <- values[[1L]]
if (!is.numeric(data$res)) cli::cli_abort("res must select a signed numeric column.")
if (is.null(title)) title <- "Process deviations"
- } else if ("osa_res" %in% names(data)) {
- if (is.null(title)) title <- "OSA residuals"
- data$res <- data$osa_res
- } else {
- if (is.null(title)) title <- "Standardized residuals"
- data$res <- data$std_res
}
data$sign <- ifelse(data$res > 0, "+", "-")
# Break size ties deterministically without advancing the caller's RNG.
@@ -439,6 +444,7 @@ plot_bubbles <- function(data, res = NULL, y = ~age, ylab = "Age", title = NULL,
if (is.null(args$colors)) {
args$colors <- c(utils::head(tam_cols(), 1), tail(tam_cols(), 1))
}
+ if (is.null(args$text)) args$text <- ~.bubble_text
p <- do.call(
plotly::plot_ly,
@@ -453,8 +459,7 @@ plot_bubbles <- function(data, res = NULL, y = ~age, ylab = "Age", title = NULL,
sizeref = 0.01,
sizemode = "area"
),
- color = ~sign,
- text = ~.bubble_text
+ color = ~sign
),
args
)
@@ -476,23 +481,11 @@ plot_resid <- function(
data,
x = ~year,
xlab = "Year",
+ residual_type = c("auto", "standardized", "osa"),
...
) {
- if ("osa_res" %in% names(data)) {
- title <- "OSA residuals"
- data$res <- data$osa_res
- } else {
- title <- "Standardized residuals"
- data$res <- data$std_res
- }
- data <- data[!is.na(data$res), ]
-
- plot_ly(data, x = x, y = ~res, ...) |>
- add_markers(alpha = 0.8) |>
- layout(
- xaxis = list(title = xlab),
- yaxis = list(title = title)
- )
+ .plot_observation_residuals(data, match.arg(residual_type), "scatter",
+ list(...), x = x, xlab = xlab)
}
#' @rdname plot_tam
diff --git a/R/population-tables.R b/R/population-tables.R
new file mode 100644
index 00000000..7c3ea870
--- /dev/null
+++ b/R/population-tables.R
@@ -0,0 +1,59 @@
+# Population-only collection shared by tidy_tam() and the quantity getters.
+# Keep comparison conversions here so lightweight consumers use the same units
+# and reported-versus-log SE conventions as dashboards.
+.population_tables <- function(model_list, interval) {
+ pop_list <- lapply(model_list, function(fit) {
+ if (inherits(fit, "tam_ref")) return(fit$pop)
+ stored <- attr(fit$pop, "interval", exact = TRUE)
+ if (!is.null(fit$pop) && !is.null(stored) &&
+ isTRUE(all.equal(stored, interval))) {
+ fit$pop
+ } else {
+ tidy_pop(fit, interval = interval)
+ }
+ })
+ references <- Filter(function(fit) inherits(fit, "tam_ref"), model_list)
+ comparison_scales <- lapply(references, `[[`, "comparison_scales")
+ comparison_scales <- Filter(function(x) length(x) > 0L, comparison_scales)
+ if (length(comparison_scales)) {
+ comparison_scales <- lapply(comparison_scales, function(x) {
+ if (is.list(x)) x <- unlist(x, use.names = TRUE)
+ if (!is.numeric(x) || is.null(names(x)) || any(!nzchar(names(x))) ||
+ anyDuplicated(names(x)) || any(!is.finite(x)) || any(x <= 0)) {
+ cli::cli_abort("A {.cls tam_ref} {.field comparison_scales} value must be a named vector of positive finite numbers.")
+ }
+ x
+ })
+ if (length(comparison_scales) > 1L &&
+ !all(vapply(comparison_scales[-1L], identical, logical(1), comparison_scales[[1L]]))) {
+ cli::cli_abort("Assessment references use different {.field comparison_scales}; the fitted outputs cannot be shown on one scale.")
+ }
+ comparison_scales <- comparison_scales[[1L]]
+ reference_pop <- references[[1L]]$pop
+ for (i in seq_along(model_list)) {
+ if (inherits(model_list[[i]], "tam_ref")) next
+ for (metric in intersect(names(comparison_scales), names(pop_list[[i]]))) {
+ tab <- pop_list[[i]][[metric]]
+ scale <- comparison_scales[[metric]]
+ reference <- reference_pop[[metric]]
+ from <- unique(tab$unit[!is.na(tab$unit)])
+ to <- unique(reference$unit[!is.na(reference$unit)])
+ if (length(from) == 1L && length(to) == 1L) .check_unit_conversion(from, to, scale)
+ for (field in intersect(c("est", "lwr", "upr"), names(tab))) {
+ tab[[field]] <- tab[[field]] * scale
+ }
+ if ("se" %in% names(tab)) {
+ natural <- if ("se_scale" %in% names(tab)) tab$se_scale == "reported" else rep(TRUE, nrow(tab))
+ tab$se[which(natural)] <- tab$se[which(natural)] * scale
+ }
+ if (is.data.frame(reference) && "unit" %in% names(reference)) {
+ units <- unique(as.character(reference$unit[!is.na(reference$unit) &
+ nzchar(as.character(reference$unit))]))
+ if (length(units) == 1L) tab$unit <- units[[1L]]
+ }
+ pop_list[[i]][[metric]] <- tab
+ }
+ }
+ }
+ pop_list
+}
diff --git a/R/prep.R b/R/prep.R
index 7f476a5d..0bb32bfe 100644
--- a/R/prep.R
+++ b/R/prep.R
@@ -662,7 +662,7 @@ prepare_tam <- function(
))
}
- check_tam_data(obs)
+ check_tam_data(obs, quiet = TRUE)
for (nm in c("years", "ages")) {
x <- get(nm)
diff --git a/R/residual-selection.R b/R/residual-selection.R
new file mode 100644
index 00000000..01556407
--- /dev/null
+++ b/R/residual-selection.R
@@ -0,0 +1,206 @@
+# Observation plots consume retained residuals; selecting a view never invokes
+# residuals() or oneStepPredict(). Missing OSA values stay missing within a model.
+.observation_residual_data <- function(data, type, selection = NULL) {
+ values <- function(name) if (name %in% names(data)) data[[name]] else rep(NA_real_, nrow(data))
+ standardized <- values("std_res")
+ osa <- values("osa_res")
+ eligible <- rep(TRUE, nrow(data))
+ if ("obs" %in% names(data)) eligible <- eligible & is.finite(data$obs) & data$obs > 0
+ if ("is_proj" %in% names(data)) eligible <- eligible & !is.na(data$is_proj) & !data$is_proj
+ osa[!eligible] <- standardized[!eligible] <- NA_real_
+ model <- if ("model" %in% names(data)) as.character(data$model) else rep("Model", nrow(data))
+ models <- if (is.factor(data$model)) levels(data$model) else unique(model)
+ models <- models[models %in% model & !is.na(models)]
+ if (!length(models)) models <- "Model"
+ method <- if ("osa_method" %in% names(data)) as.character(data$osa_method) else rep(NA_character_, nrow(data))
+ method[is.na(method) | !nzchar(method)] <- "method unspecified"
+ data$.residual_value <- rep(NA_real_, nrow(data))
+ data$.residual_type <- rep(NA_character_, nrow(data))
+ data$.residual_source <- rep(NA_character_, nrow(data))
+ status <- lapply(models, function(name) {
+ rows <- which(model == name)
+ selected <- if (type == "auto" && !is.null(selection)) selection[[name]] else if (type == "auto") {
+ if (any(is.finite(osa[rows]))) "osa" else "standardized"
+ } else type
+ data$.residual_value[rows] <<- if (selected == "osa") osa[rows] else standardized[rows]
+ data$.residual_type[rows] <<- selected
+ label <- if (selected == "osa") paste0("OSA residuals (", method[rows], ")") else "Standardized residuals"
+ data$.residual_source[rows] <<- label
+ finite <- rows[is.finite(data$.residual_value[rows])]
+ data.frame(model = name, residual_type = selected,
+ osa_method = if (selected == "osa" && length(finite)) paste(unique(method[finite]), collapse = ", ") else NA_character_,
+ n = length(finite), stringsAsFactors = FALSE)
+ })
+ list(data = data, status = do.call(rbind, status))
+}
+
+.residual_axis_range <- function(data, formula) {
+ value <- stats::model.frame(formula, data, na.action = stats::na.pass)[[1L]]
+ if (!is.numeric(value)) return(NULL)
+ value <- value[is.finite(value)]
+ if (!length(value)) return(NULL)
+ limits <- range(value)
+ limits + c(-1, 1) * max(diff(limits) * .05, .1)
+}
+
+.residual_frame_values <- function(data, frame) {
+ if (is.null(frame)) return(list(values = rep("Model", nrow(data)), levels = "Model", prefix = ""))
+ value <- if (inherits(frame, "formula")) {
+ tab <- stats::model.frame(frame, data, na.action = stats::na.pass)
+ if (ncol(tab) != 1L) cli::cli_abort("Residual plot frames must select one column.")
+ tab[[1L]]
+ } else if (is.character(frame) && length(frame) == 1L && frame %in% names(data)) {
+ data[[frame]]
+ } else frame
+ if (length(value) != nrow(data) || anyNA(value)) {
+ cli::cli_abort("Residual plot frames must provide one nonmissing value per row.")
+ }
+ levels <- if (is.factor(value)) levels(value)[levels(value) %in% value] else unique(as.character(value))
+ if (!length(levels)) levels <- "Model"
+ label <- if (inherits(frame, "formula")) all.vars(frame)[1L] else "Frame"
+ if (!length(label) || is.na(label)) label <- "Frame"
+ list(values = as.character(value), levels = levels,
+ prefix = paste0(toupper(substr(label, 1L, 1L)), substring(label, 2L), ": "))
+}
+
+.residual_plot_label <- function(status, type) {
+ available <- status[status$n > 0L, , drop = FALSE]
+ kinds <- unique(available$residual_type)
+ if (identical(kinds, "osa")) "OSA residuals" else if (identical(kinds, "standardized")) {
+ "Standardized residuals"
+ } else if (type == "osa") "OSA residuals" else if (type == "standardized") {
+ "Standardized residuals"
+ } else "Residuals (Auto)"
+}
+
+.residual_plot_annotations <- function(status, type) {
+ available <- status$n > 0L
+ if (!any(available)) {
+ text <- if (type == "osa" || all(status$residual_type == "osa")) "OSA residuals unavailable for this model/panel" else
+ "Residuals unavailable for this model/panel"
+ return(list(list(text = text, x = .5, y = .5, xref = "paper", yref = "paper", showarrow = FALSE)))
+ }
+ labels <- ifelse(status$residual_type == "osa", paste0("OSA (", status$osa_method, ")"), "Standardized")
+ labels[!available] <- if (type == "osa") "OSA unavailable" else "Unavailable"
+ if (nrow(status) > 1L) labels <- paste0(status$model, ": ", labels)
+ text <- paste(labels, collapse = "; ")
+ if (type == "auto") text <- paste0("Auto: ", text)
+ list(list(text = as.character(htmltools::htmlEscape(text)), x = 0, y = 1,
+ xref = "paper", yref = "paper", xanchor = "left", yanchor = "bottom",
+ yshift = 8, showarrow = FALSE, font = list(size = 10)))
+}
+
+.residual_plot_cell <- function(data, type, kind, args, x, xlab, y, ylab, title, selection) {
+ selected <- .observation_residual_data(data, type, selection)
+ d <- selected$data
+ label <- .residual_plot_label(selected$status, type)
+ d <- d[is.finite(d$.residual_value), , drop = FALSE]
+ if (!nrow(d)) {
+ p <- plotly::plot_ly(type = "scatter", mode = "markers", x = numeric(), y = numeric())
+ } else if (kind == "bubbles") {
+ d$kind <- as.character(htmltools::htmlEscape(d$.residual_source))
+ p <- do.call(plot_bubbles, c(list(data = d, res = ~.residual_value,
+ y = y, ylab = ylab, title = if (is.null(title)) label else title), args))
+ } else {
+ d$.residual_hover <- paste0(htmltools::htmlEscape(d$.residual_source), ": ", signif(d$.residual_value, 4))
+ if (is.null(args$text)) args$text <- ~.residual_hover
+ if (is.null(args$hovertemplate)) args$hovertemplate <- "%{x}
%{text}%{fullData.name}"
+ p <- do.call(plotly::plot_ly, c(list(data = d, x = x, y = ~.residual_value), args)) |>
+ plotly::add_markers(alpha = .8)
+ }
+ p <- plotly::plotly_build(p)
+ if (!nrow(d)) for (j in seq_along(p$x$data)) {
+ p$x$data[[j]]$x <- p$x$data[[j]]$y <- numeric()
+ p$x$data[[j]]$showlegend <- FALSE
+ }
+ p$x$layout$xaxis$title <- if (kind == "bubbles") "Year" else xlab
+ p$x$layout$yaxis$title <- if (kind == "bubbles") ylab else label
+ if (kind == "bubbles") p$x$layout$title <- if (is.null(title)) label else title
+ p$x$layout$annotations <- .residual_plot_annotations(selected$status, type)
+ list(data = p$x$data, layout = p$x$layout, status = selected$status)
+}
+
+.residual_trace_key <- function(trace) {
+ paste(trace$name, trace$legendgroup, trace$type, trace$mode, collapse = "\r")
+}
+
+.plot_observation_residuals <- function(data, type, kind, args,
+ x = ~year, xlab = "Year", y = ~age,
+ ylab = "Age", title = NULL) {
+ if (!is.data.frame(data)) cli::cli_abort("Residual plots require a data frame.")
+ frame <- .residual_frame_values(data, args$frame)
+ animated <- !is.null(args$frame) && length(frame$levels) > 1L
+ args$frame <- NULL
+ types <- c("auto", "standardized", "osa")
+ models <- if (is.factor(data$model)) levels(data$model) else unique(as.character(data$model))
+ model_colors <- NULL
+ if (kind == "scatter" && inherits(args$color, "formula") && identical(all.vars(args$color), "model")) {
+ palette <- args$colors
+ if (is.null(palette)) palette <- tam_pal(length(models))
+ if (is.function(palette)) palette <- palette(length(models))
+ if (is.character(palette) && length(palette) >= length(models)) model_colors <- stats::setNames(palette[seq_along(models)], models)
+ }
+ sign_colors <- if (is.null(args$colors)) c(tam_cols()[1L], tail(tam_cols(), 1L)) else args$colors
+ if (is.function(sign_colors)) sign_colors <- sign_colors(2L)
+ views <- lapply(types, function(view) {
+ chosen <- .observation_residual_data(data, view)
+ selection <- stats::setNames(chosen$status$residual_type, chosen$status$model)
+ chosen <- chosen$data
+ valid <- chosen[is.finite(chosen$.residual_value), , drop = FALSE]
+ xrange <- .residual_axis_range(valid, if (kind == "bubbles") ~year else x)
+ yrange <- .residual_axis_range(valid, if (kind == "bubbles") y else ~.residual_value)
+ cells <- lapply(frame$levels, function(level) {
+ d <- if (animated) data[frame$values == level, , drop = FALSE] else data
+ cell <- .residual_plot_cell(d, view, kind, args, x, xlab, y, ylab, title, selection)
+ if (!is.null(xrange)) cell$layout$xaxis <- utils::modifyList(cell$layout$xaxis, list(range = xrange, autorange = FALSE))
+ if (!is.null(yrange)) cell$layout$yaxis <- utils::modifyList(cell$layout$yaxis, list(range = yrange, autorange = FALSE))
+ for (j in seq_along(cell$data)) {
+ trace <- cell$data[[j]]
+ if (kind == "bubbles" && length(trace$name) == 1L &&
+ trace$name %in% c("+", "-") && is.character(sign_colors)) {
+ trace$marker$color <- sign_colors[match(trace$name, c("+", "-"))]
+ } else if (length(model_colors) && length(trace$name) == 1L && trace$name %in% names(model_colors)) {
+ trace$marker$color <- unname(model_colors[trace$name])
+ }
+ cell$data[[j]] <- trace
+ }
+ cell
+ })
+ list(type = view, cells = cells)
+ })
+ # Stable trace slots clear markers when a model lacks the selected residuals.
+ templates <- list()
+ for (view in views) for (cell in view$cells) for (trace in cell$data) {
+ key <- .residual_trace_key(trace)
+ if (is.null(templates[[key]])) templates[[key]] <- trace
+ }
+ for (i in seq_along(views)) for (j in seq_along(views[[i]]$cells)) {
+ cell <- views[[i]]$cells[[j]]
+ keys <- vapply(cell$data, .residual_trace_key, character(1))
+ cell$data <- lapply(names(templates), function(key) {
+ found <- match(key, keys)
+ if (!is.na(found)) return(cell$data[[found]])
+ trace <- templates[[key]]
+ trace$x <- trace$y <- numeric()
+ trace$text <- character()
+ trace$showlegend <- FALSE
+ if (!is.null(trace$marker$size)) trace$marker$size <- numeric()
+ trace
+ })
+ views[[i]]$cells[[j]] <- cell
+ }
+ config <- list(views = views, frames = frame$levels, animated = animated,
+ prefix = frame$prefix, view = match(type, types) - 1L, frame = 0L)
+ initial <- views[[match(type, types)]]$cells[[1L]]
+ p <- plotly::plotly_build(plotly::plot_ly(type = "scatter", mode = "markers", x = numeric(), y = numeric()))
+ p$x$data <- initial$data
+ p$x$layout <- utils::modifyList(p$x$layout, initial$layout)
+ if (animated) {
+ cells <- views[[match(type, types)]]$cells
+ p$x$frames <- lapply(seq_along(frame$levels), function(i) list(
+ name = frame$levels[i], data = cells[[i]]$data, layout = cells[[i]]$layout,
+ traces = seq_along(cells[[i]]$data) - 1L))
+ }
+ script <- paste(readLines(system.file("js", "residual-selection.js", package = "tinyAM"), warn = FALSE), collapse = "\n")
+ htmlwidgets::onRender(p, script, data = config)
+}
diff --git a/R/residuals.R b/R/residuals.R
index 967bd82a..4a6d728a 100644
--- a/R/residuals.R
+++ b/R/residuals.R
@@ -6,50 +6,52 @@
#' @param object A fitted `tam_fit`.
#' @param type `"standardized"` (default), `"log"` (log observation minus log
#' prediction), `"response"` (observation minus prediction), or `"osa"`.
-#' @param ... Computation controls passed to [tidy_obs_pred()] and
-#' [RTMB::oneStepPredict()] when uncached OSA residuals are requested, such as
+#' @param method OSA algorithm, as in [add_osa_residuals()]. Only used with
+#' `type = "osa"`. Explicitly supplying it requests a fresh calculation.
+#' @param ... Computation controls passed to [RTMB::oneStepPredict()], such as
#' `trace = FALSE` or `parallel = FALSE`.
#' `subset` and `conditional` index included observations in likelihood order;
#' unrequested rows retain missing residuals.
-#' Controls are not used for cached OSA values; other residual types reject them.
-#' @details `type = "osa"` reuses cached OSA columns, or explicitly requests
-#' their calculation at fitted parameter values when absent. It can be expensive,
+#' Explicit controls request a fresh calculation; other residual types reject them.
+#' @details `type = "osa"` reuses cached OSA columns when no computation controls
+#' or method are supplied. Otherwise both catch and index residuals are calculated
+#' together, using one algorithm at fitted parameter values. It can be expensive,
#' but does not rerun [fit_tam()] or re-estimate fixed parameters. Temporary
#' objectives protect the fitted object's state. Other residual types remain
#' simple extractors. Excluded observations, including zeros, retain rows with
#' missing residuals. The supplied fit is not modified or implicitly cached.
-#' To retain OSA results for later extraction and dashboard use, assign
-#' `fit$obs_pred <- tidy_obs_pred(fit, add_osa_res = TRUE, trace = FALSE)`.
+#' To retain results for later extraction and dashboard use, assign
+#' `fit_osa <- add_osa_residuals(fit, method = "oneStepGeneric", trace = FALSE)`.
#' @return A data frame with year, age, type (catch/index), survey, observations,
-#' predictions, residuals, and available unit/fleet columns. The `residual_type`
-#' attribute names the requested scale.
+#' predictions, residuals, and available unit/fleet columns. OSA output also
+#' identifies `osa_method`. The `residual_type` attribute names the requested
+#' scale, and `osa` retains available method/settings/coverage metadata.
#' @examples
#' \dontrun{
#' residuals(fit)
#' residuals(fit, type = "log")
-#' residuals(fit, type = "osa", trace = FALSE)
+#' residuals(fit, type = "osa", method = "oneStepGeneric", trace = FALSE)
#' }
#' @seealso [tidy_obs_pred()], [plot_bubbles()]
#' @importFrom stats residuals
#' @export
-residuals.tam_fit <- function(object, type = c("standardized", "log", "response", "osa"), ...) {
+residuals.tam_fit <- function(object, type = c("standardized", "log", "response", "osa"),
+ method = "oneStepGaussianOffMode", ...) {
object <- .require_tam_fit(object, arg = "object")
type <- match.arg(type)
- tables <- object$obs_pred[c("catch", "index")]
- if (type != "osa" && length(list(...))) cli::cli_abort("Additional computation controls are only used with type = 'osa'.")
+ controls <- list(...)
+ explicit <- !missing(method) || length(controls) > 0L
+ .osa_method(method, controls)
+ if (type != "osa" && explicit) cli::cli_abort("Additional computation controls are only used with type = 'osa'.")
+ tables <- .observation_tables(object)
+ metadata <- attr(tables, "osa", exact = TRUE)
+ if (is.null(metadata)) metadata <- object$osa
missing_osa <- vapply(tables, function(d) is.data.frame(d) && nrow(d) > 0L &&
!"osa_res" %in% names(d) && any(is.finite(d$obs) & d$obs > 0 & if (is.null(d$is_proj)) TRUE else !d$is_proj), logical(1))
- if (type == "osa" && any(missing_osa)) {
- computed <- tryCatch(tidy_obs_pred(object, add_osa_res = TRUE, ...),
- error = function(e) cli::cli_abort("Could not calculate OSA residuals: {conditionMessage(e)}", parent = e))
- for (name in names(tables)[missing_osa]) {
- d <- tables[[name]]; full <- computed[[name]]
- keys <- intersect(c("year", "age", "survey", "fleet"), names(d))
- key <- function(x) do.call(paste, c(x[keys], sep = "\r"))
- i <- match(key(d), key(full))
- if (anyNA(i)) cli::cli_abort("OSA observation rows do not match the fitted inputs.")
- tables[[name]]$osa_res <- full$osa_res[i]
- }
+ if (type == "osa" && (explicit || any(missing_osa))) {
+ computed <- .compute_osa_residuals(object, method, controls)
+ metadata <- computed$metadata
+ tables <- .attach_osa_tables(tables, computed$tables, metadata)
}
out <- .bind_rows(lapply(names(tables), function(name) {
@@ -68,12 +70,18 @@ residuals.tam_fit <- function(object, type = c("standardized", "log", "response"
keep <- is.finite(d$obs) & d$obs > 0 & is.finite(residual)
if ("is_proj" %in% names(d)) keep <- keep & !d$is_proj
residual[!keep] <- NA_real_
- d <- d[intersect(c("year", "age", "survey", "fleet", "obs", "pred", "unit", "is_proj"), names(d))]
+ if (type == "osa" && !"osa_method" %in% names(d)) {
+ d$osa_method <- rep(if (is.null(metadata$method)) NA_character_ else metadata$method, nrow(d))
+ }
+ columns <- c("year", "age", "survey", "fleet", "obs", "pred", "unit", "is_proj")
+ if (type == "osa") columns <- c(columns, "osa_method")
+ d <- d[intersect(columns, names(d))]
d$type <- name; d$residual <- residual
if (!"survey" %in% names(d)) d$survey <- NA_character_
d
}))
attr(out, "residual_type") <- type
+ if (type == "osa" && !is.null(metadata)) attr(out, "osa") <- metadata
out
}
diff --git a/R/smooths.R b/R/smooths.R
index 916c036f..98b733a5 100644
--- a/R/smooths.R
+++ b/R/smooths.R
@@ -175,9 +175,90 @@ NULL
list(fixed=fixed,terms=terms)
}
+# A smooth2random basis can change eigenvector signs (or rotate an eigenspace)
+# without changing its function or Gaussian penalty. Transfer coefficients only
+# when both the function and, for random coefficients, the penalty are preserved.
+.smooth_design_transform <- function(old, new, orthogonal = TRUE) {
+ tolerance <- 1e-10
+ if (!is.matrix(old) || !is.matrix(new) || !identical(dim(old), dim(new)) ||
+ any(!is.finite(old)) || any(!is.finite(new))) return(NULL)
+ n <- ncol(old)
+ if (!n) return(matrix(numeric(), 0L, 0L))
+ if (identical(unname(old), unname(new))) return(diag(n))
+ transform <- tryCatch(qr.solve(new, old, tol = tolerance),
+ error = function(e) NULL)
+ if (is.null(transform) || any(!is.finite(transform))) return(NULL)
+ scale <- max(abs(old))
+ if (!scale || max(abs(new %*% transform - old)) > tolerance * scale) return(NULL)
+ if (orthogonal && max(abs(crossprod(transform) - diag(n))) > tolerance) return(NULL)
+ if (!orthogonal && qr(transform, tol = tolerance)$rank != n) return(NULL)
+ transform
+}
+
+.smooth_basis_transform <- function(old, new) {
+ fields <- c("type", "smooth_id", "variables", "by", "scope")
+ if (!identical(old$type, "smooth") || !identical(new$type, "smooth") ||
+ !identical(old[fields], new[fields]) ||
+ !identical(class(old$smooth), class(new$smooth)) ||
+ !identical(old$smooth$bs.dim, new$smooth$bs.dim) ||
+ !isTRUE(all.equal(old$basis_signature$grid, new$basis_signature$grid,
+ tolerance = 0))) return(NULL)
+ before <- old$basis_signature$design
+ after <- new$basis_signature$design
+ if (length(before$random) != length(after$random)) return(NULL)
+ fixed <- .smooth_design_transform(before$fixed, after$fixed, orthogonal = FALSE)
+ random <- Map(.smooth_design_transform, before$random, after$random)
+ if (is.null(fixed) || any(vapply(random, is.null, logical(1)))) return(NULL)
+ list(fixed = fixed, random = random)
+}
+
.smooth_basis_compatible <- function(old, new) {
- identical(old$type,"smooth") && identical(new$type,"smooth") &&
- isTRUE(all.equal(old$basis_signature,new$basis_signature,tolerance=1e-10))
+ !is.null(.smooth_basis_transform(old, new))
+}
+
+.smooth_warm_starts <- function(start, old, new) {
+ designs <- list(F = old$F_modmat, q = old$q_modmat)
+ for (component in c("q", "F")) {
+ field <- if (component == "q") "q_terms" else "F_terms"
+ coefficient <- if (component == "F") "log_mu_f" else
+ if (new$index_settings$q_link == "logit") "logit_q" else "log_q"
+ terms <- Filter(function(x) x$type == "smooth", new[[field]])
+ ids <- unique(vapply(terms, `[[`, character(1), "smooth_id"))
+ for (id in ids) {
+ pieces <- Filter(function(x) identical(x$smooth_id, id), terms)
+ previous <- Filter(function(x) identical(x$type, "smooth") &&
+ identical(x$smooth_id, id), old[[field]])
+ transform <- if (identical(names(previous), names(pieces)))
+ .smooth_basis_transform(previous[[1L]], pieces[[1L]]) else NULL
+ columns <- colnames(pieces[[1L]]$smooth_fixed)
+ if (is.null(transform)) {
+ parameters <- unlist(lapply(pieces, function(x) c(x$parameter, x$sd_parameter)))
+ start[parameters] <- NULL
+ if (!is.null(start[[coefficient]])) {
+ start[[coefficient]] <- start[[coefficient]][!names(start[[coefficient]]) %in% columns]
+ }
+ next
+ }
+ for (piece in pieces) {
+ parameter <- piece$parameter
+ if (!is.null(start[[parameter]])) {
+ value <- drop(transform$random[[piece$penalty]] %*% start[[parameter]])
+ names(value) <- names(start[[parameter]])
+ start[[parameter]] <- value
+ }
+ }
+ if (length(columns)) {
+ if (!is.null(start[[coefficient]]) && all(columns %in% names(start[[coefficient]]))) {
+ start[[coefficient]][columns] <- drop(transform$fixed %*% start[[coefficient]][columns])
+ }
+ # Express the old design in the new coordinates before checking the
+ # ordinary fixed-effect design, including the smooth null space.
+ designs[[component]][, columns] <- designs[[component]][, columns, drop = FALSE] %*%
+ solve(transform$fixed)
+ }
+ }
+ }
+ list(parameters = start, designs = designs)
}
.smooth_function_values <- function(par, dat, component, effects) {
diff --git a/R/tidy.R b/R/tidy.R
index e856501a..ab4a0b97 100644
--- a/R/tidy.R
+++ b/R/tidy.R
@@ -73,9 +73,9 @@ tidy_mat <- tidy_array
#'
#' @param fit A fitted TAM object as returned by [fit_tam()].
#' @param interval Confidence level for combined catchability intervals.
-#' @param add_osa_res Logical; add one-step-ahead residuals? Hard wired to
-#' apply the `"oneStepGaussianOffMode"` method.
-#' See [RTMB::oneStepPredict()] for details.
+#' @param add_osa_res Logical; calculate one-step-ahead residuals?
+#' @param method OSA algorithm; defaults to `"oneStepGaussianOffMode"`.
+#' See [add_osa_residuals()] for supported choices and computation details.
#' @param ... Arguments to pass to [RTMB::oneStepPredict()].
#'
#' @details
@@ -87,6 +87,10 @@ tidy_mat <- tidy_array
#' different predictive calculation. See [tinyAM-model] for the observation equations.
#' OSA calculation rebuilds a temporary objective at stored fitted parameters;
#' it does not refit fixed parameters or modify the fitted objective.
+#' With `add_osa_res = TRUE`, it computes the requested method afresh and records
+#' its settings and observation coverage. `osa_method` identifies the algorithm
+#' beside `osa_res`. To retain diagnostics while preserving existing observation
+#' tables and confidence intervals, use [add_osa_residuals()].
#' When stockassessment is loaded, one-step residuals run sequentially: the
#' upstream TMB parallel helper cannot select between SAM and RTMB libraries.
#'
@@ -112,9 +116,10 @@ tidy_mat <- tidy_array
#' @seealso [fit_tam()], [tidy_rep()], [tidy_sdrep()], [tidy_pop()]
#' @export
-tidy_obs_pred <- function(fit, add_osa_res = FALSE, interval = .95, ...) {
+tidy_obs_pred <- function(fit, add_osa_res = FALSE, interval = .95,
+ method = "oneStepGaussianOffMode", ...) {
fit <- .require_tam_fit(fit, arg = "fit")
- osa_obj <- if (add_osa_res) .fitted_objective(fit) else NULL
+ .osa_method(method, list(...))
if (!is.numeric(interval) || length(interval) != 1L || !is.finite(interval) ||
interval <= 0 || interval >= 1) {
cli::cli_abort("{.arg interval} must be one number strictly between zero and one.")
@@ -169,26 +174,8 @@ tidy_obs_pred <- function(fit, add_osa_res = FALSE, interval = .95, ...) {
obs_pred$index$std_res <- with(obs_pred$index, ifelse(obs == 0, NA, (log(obs) - log_pred$index) / sd))
if (add_osa_res) {
- args <- list(...)
- if (isTRUE(args$parallel) && "stockassessment" %in% loadedNamespaces()) {
- args$parallel <- FALSE
- cli::cli_inform("Calculating one-step residuals sequentially while SAM and RTMB are loaded.")
- }
- osa_res <- do.call(
- RTMB::oneStepPredict,
- c(list(obj = osa_obj, method = "oneStepGaussianOffMode"), args)
- )
- selected <- if (is.null(args$subset)) setdiff(seq_along(fit$dat$observed), args$conditional) else args$subset
- if (length(osa_res$residual) != length(selected)) {
- cli::cli_abort("OSA results do not match the requested observation subset.")
- }
- all_osa <- rep(NA_real_, length(fit$dat$observed))
- all_osa[selected] <- osa_res$residual
- split_osa_res <- split(all_osa, fit$dat$obs_map$type[fit$dat$obs_map$is_observed])
- split_is_observed <- split(fit$dat$obs_map$is_observed, fit$dat$obs_map$type)
- obs_pred$catch$osa_res <- obs_pred$index$osa_res <- NA
- obs_pred$catch$osa_res[split_is_observed$catch] <- split_osa_res$catch
- obs_pred$index$osa_res[split_is_observed$index] <- split_osa_res$index
+ osa <- .compute_osa_residuals(fit, method, list(...))
+ obs_pred <- .attach_osa_tables(obs_pred, osa$tables, osa$metadata)
}
attr(obs_pred, "interval") <- interval
@@ -1048,66 +1035,12 @@ tidy_tam <- function(..., model_list = NULL, interval = 0.95, label = "model", l
return(fit$obs_pred)
}
tables <- tidy_obs_pred(fit, interval = interval)
- for (nm in names(tables)) {
- extras <- setdiff(names(fit$obs_pred[[nm]]), names(tables[[nm]]))
- tables[[nm]][extras] <- fit$obs_pred[[nm]][extras]
- }
- tables
+ metadata <- attr(fit$obs_pred, "osa", exact = TRUE)
+ if (is.null(metadata)) metadata <- fit$osa
+ .preserve_observation_extras(tables, fit$obs_pred, metadata)
})
- pop_list <- lapply(model_list, function(fit) {
- if (inherits(fit, "tam_ref")) {
- return(fit$pop)
- }
- if (!is.null(fit$pop) && interval_matches(fit$pop)) {
- fit$pop
- } else {
- tidy_pop(fit, interval = interval)
- }
- })
- references <- Filter(function(fit) inherits(fit, "tam_ref"), model_list)
- comparison_scales <- lapply(references, `[[`, "comparison_scales")
- comparison_scales <- Filter(function(x) length(x) > 0L, comparison_scales)
- if (length(comparison_scales)) {
- comparison_scales <- lapply(comparison_scales, function(x) {
- if (is.list(x)) x <- unlist(x, use.names = TRUE)
- if (!is.numeric(x) || is.null(names(x)) || any(!nzchar(names(x))) ||
- anyDuplicated(names(x)) || any(!is.finite(x)) || any(x <= 0)) {
- cli::cli_abort("A {.cls tam_ref} {.field comparison_scales} value must be a named vector of positive finite numbers.")
- }
- x
- })
- if (length(comparison_scales) > 1L &&
- !all(vapply(comparison_scales[-1L], identical, logical(1), comparison_scales[[1L]]))) {
- cli::cli_abort("Assessment references use different {.field comparison_scales}; the fitted outputs cannot be shown on one scale.")
- }
- comparison_scales <- comparison_scales[[1L]]
- reference_pop <- references[[1L]]$pop
- for (i in seq_along(model_list)) {
- if (inherits(model_list[[i]], "tam_ref")) next
- for (metric in intersect(names(comparison_scales), names(pop_list[[i]]))) {
- tab <- pop_list[[i]][[metric]]
- scale <- comparison_scales[[metric]]
- reference <- reference_pop[[metric]]
- from <- unique(tab$unit[!is.na(tab$unit)])
- to <- unique(reference$unit[!is.na(reference$unit)])
- if (length(from) == 1L && length(to) == 1L) .check_unit_conversion(from, to, scale)
- for (field in intersect(c("est", "lwr", "upr"), names(tab))) {
- tab[[field]] <- tab[[field]] * scale
- }
- if ("se" %in% names(tab)) {
- natural <- if ("se_scale" %in% names(tab)) tab$se_scale == "reported" else rep(TRUE, nrow(tab))
- tab$se[which(natural)] <- tab$se[which(natural)] * scale
- }
- if (is.data.frame(reference) && "unit" %in% names(reference)) {
- units <- unique(as.character(reference$unit[!is.na(reference$unit) &
- nzchar(as.character(reference$unit))]))
- if (length(units) == 1L) tab$unit <- units[[1L]]
- }
- pop_list[[i]][[metric]] <- tab
- }
- }
- }
+ pop_list <- .population_tables(model_list, interval)
fitted <- which(vapply(model_list, inherits, logical(1), "tam_fit"))
pattern_units <- if (length(fitted)) {
displayed <- pop_list[[fitted[1L]]]
@@ -1135,6 +1068,17 @@ tidy_tam <- function(..., model_list = NULL, interval = 0.95, label = "model", l
})
obs_pred <- stack_nested(obs_list, label = id_col, label_type = label_type)
+ osa_metadata <- lapply(seq_along(model_list), function(i) {
+ metadata <- attr(obs_list[[i]], "osa", exact = TRUE)
+ if (is.null(metadata)) metadata <- model_list[[i]]$osa
+ metadata
+ })
+ names(osa_metadata) <- names(model_list)
+ osa_metadata <- Filter(Negate(is.null), osa_metadata)
+ if (length(osa_metadata)) {
+ attr(obs_pred, "osa") <- osa_metadata
+ for (name in names(obs_pred)) attr(obs_pred[[name]], "osa") <- osa_metadata
+ }
pop <- stack_nested(pop_list, label = id_col, label_type = label_type)
.unit_label(obs_pred$catch, "")
if (is.data.frame(obs_pred$index)) {
@@ -1154,6 +1098,7 @@ tidy_tam <- function(..., model_list = NULL, interval = 0.95, label = "model", l
out <- list(
obs_pred = obs_pred,
+ osa = osa_metadata,
pop = pop,
patterns = stack_nested(patterns, label = id_col, label_type = label_type),
fixed_par = fixed_tbl,
diff --git a/R/vis.R b/R/vis.R
index 11c00a43..e1c45e5b 100644
--- a/R/vis.R
+++ b/R/vis.R
@@ -21,6 +21,9 @@
#' form a uniquely named list so they can be labeled in the dashboard.
#' Additional arguments for [rmarkdown::render()] can be passed through
#' `render_args`, which must itself be a (named) list.
+#' Dashboard R packages are installed with tinyAM. Rendering also needs
+#' Pandoc, included with RStudio or available from
+#' .
#' Precomputed assessment references may be mixed with tinyAM fits; they are
#' not refittable TAM objects. Supply `background` to add a page
#' describing assessment assumptions and translation choices.
@@ -31,6 +34,11 @@
#' Without fleets, Aggregate fishery retains all single-fishery
#' diagnostics. With fleets, it shows total yield; residuals are
#' assessed separately within each fleet.
+#' Observation residuals default to Auto: each model uses retained OSA
+#' residuals when available in that panel, otherwise standardized
+#' residuals. The dropdown also selects either type explicitly. Partial
+#' OSA gaps remain missing, and labels identify the method. Retain OSA
+#' calculations with [add_osa_residuals()]; rendering never computes them.
#' The Parameters menu separates Fixed and Random pages. Fixed plots
#' group similar quantities, with initial abundance shown separately.
#' Stock-recruit models add a relationship page with parent-year pairs,
@@ -77,6 +85,8 @@
#' }
#'
#' @importFrom utils browseURL
+#' @importFrom knitr knit_expand
+#' @importFrom flexdashboard flex_dashboard
#'
#' @export
vis_tam <- function(..., model_list = NULL, interval = 0.95, output_file = NULL,
@@ -190,6 +200,12 @@ vis_tam.tam_retro <- function(..., model_list = NULL, interval = 0.95, output_fi
"i" = "Install with: {.code {install_call}}"
))
}
+ if (!rmarkdown::pandoc_available()) {
+ cli::cli_abort(c(
+ "Pandoc is required to render a dashboard.",
+ "i" = "RStudio includes Pandoc. Otherwise install it from https://pandoc.org/installing.html and make it available on PATH."
+ ))
+ }
if (!is.null(background)) {
if (!is.character(background) || !length(background) || anyNA(background)) {
diff --git a/README.md b/README.md
index 005ca102..fb47fc10 100644
--- a/README.md
+++ b/README.md
@@ -14,9 +14,12 @@ The development version can be installed from GitHub:
```r
# install.packages("pak")
-pak::pak(c("PaulRegular/tinyAM", "rmarkdown", "flexdashboard"))
+pak::pak("PaulRegular/tinyAM")
```
+Dashboard R packages install automatically. Rendering requires
+[Pandoc](https://pandoc.org/installing.html), which is included with RStudio.
+
See the [package website](https://paulregular.github.io/tinyAM/) and
[cheatsheet](https://paulregular.github.io/tinyAM/cheatsheet/tinyAM_cheatsheet.pdf)
for guidance and examples.
@@ -212,8 +215,15 @@ statistical advisories. Convergence alone does not establish that the data
identify every component.
Inspect each check before interpreting or ranking an alternative.
+Use `plot_model(base, N_iid, M_ar1)` or the dashboard's **Model structure** page
+to inspect process switches and formulas. The
+[structure guide](https://paulregular.github.io/tinyAM/articles/model-structure.html)
+links contrasting diagrams to the model equations.
+
In the dashboard, start with **Fisheries** and **Surveys** to inspect fits and
-standardized residuals. Fisheries lists each fleet, then **Aggregate fishery**;
+residuals. **Auto** uses retained OSA residuals where available for each model;
+otherwise it uses standardized residuals. Select either type explicitly with
+the residual dropdown. Fisheries lists each fleet, then **Aggregate fishery**;
without fleets, the aggregate page contains the single-fishery diagnostics.
Review **Tables → Diagnostics** and **Parameters** for
checks, uncertainty and process deviations before interpreting **Trends** for
@@ -251,8 +261,12 @@ before interpreting rankings.
`tidy_tam(base, N_iid, M_ar1)` collects population, parameter and observation tables.
Confidence limits use the reported scale; `se_scale` identifies the SE scale.
- [Residuals](https://paulregular.github.io/tinyAM/reference/residuals.tam_fit.html):
- `residuals(base)` extracts standardized residuals; `type = "osa"` explicitly
- requests one-step-ahead calculations without refitting.
+ `residuals(base)` extracts standardized residuals. To calculate and retain
+ OSAs for dashboards, use `base <- add_osa_residuals(base, method = "oneStepGeneric")`.
+ Calculation can be expensive and does not refit the model.
+- [Population shortcuts](https://paulregular.github.io/tinyAM/reference/get_population.html):
+ `get_ssb(base)`, `get_recruitment(base)`, `get_fbar(base)` and `get_mbar(base)`
+ select annual summaries; `get_f(base)` and `get_m(base)` select age-specific tables.
- [Simulation](https://paulregular.github.io/tinyAM/reference/sim_tam.html) and
[projections](https://paulregular.github.io/tinyAM/reference/prepare_tam.html)
explore model behavior and future scenarios.
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 708c30aa..9b741c35 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -13,12 +13,15 @@ template:
heading_font: 'Arial, sans-serif'
navbar:
structure:
- left: [intro, reference, cheatsheet]
+ left: [intro, model_structure, reference, cheatsheet]
right: [search, github]
components:
intro:
text: Home
href: index.html
+ model_structure:
+ text: Model structure
+ href: articles/model-structure.html
cheatsheet:
text: Cheatsheet
href: cheatsheet/tinyAM_cheatsheet.pdf
@@ -60,6 +63,7 @@ reference:
- title: Visualize and compare
- contents:
- vis_tam
+ - plot_model
- plot_tam
- plot_age_year
- plot_matrix
@@ -70,7 +74,9 @@ reference:
- tidy_tam
- tidy_obs_pred
- residuals.tam_fit
+ - add_osa_residuals
- tidy_pop
+ - get_population
- tidy_par
- tidy_recruitment
- tidy_sdrep
diff --git a/inst/js/residual-selection.js b/inst/js/residual-selection.js
new file mode 100644
index 00000000..8bed14f1
--- /dev/null
+++ b/inst/js/residual-selection.js
@@ -0,0 +1,63 @@
+function(el, x, config) {
+ 'use strict';
+ function array(value) { return Array.isArray(value) ? value : [value]; }
+ config.frames = array(config.frames);
+ config.views.forEach(function(view) {
+ view.cells.forEach(function(cell) {
+ cell.data.forEach(function(trace) {
+ ['x', 'y', 'text', 'hovertext', 'customdata'].forEach(function(key) {
+ if (trace[key] !== undefined) trace[key] = array(trace[key]);
+ });
+ if (trace.marker && trace.marker.size !== undefined) trace.marker.size = array(trace.marker.size);
+ });
+ });
+ });
+ var state = {view: config.view, frame: config.frame}, busy = false, pending = false;
+ function draw() {
+ if (busy) { pending = true; return; }
+ busy = true;
+ var view = config.views[state.view], cell = view.cells[state.frame];
+ var frames = config.frames.map(function(name, i) {
+ return {name: name, data: view.cells[i].data, layout: view.cells[i].layout,
+ traces: view.cells[i].data.map(function(trace, j) { return j; })};
+ });
+ var layout = Object.assign({}, x.layout, cell.layout, {
+ autosize: true, margin: Object.assign({}, cell.layout.margin, {t: 100}),
+ // Leave the upper-right modebar clear and reserve a separate method row.
+ updatemenus: [{name: 'tam-residual-type', type: 'dropdown',
+ x: 0, y: 1, xanchor: 'left', yanchor: 'bottom', pad: {b: 28}, active: state.view,
+ buttons: ['Auto', 'Standardized', 'OSA'].map(function(label, i) {
+ return {label: label, method: 'skip', execute: false, args: [i]};
+ })}],
+ sliders: config.animated ? [{active: state.frame, pad: {t: 45},
+ currentvalue: {prefix: config.prefix}, steps: config.frames.map(function(name) {
+ return {label: name, method: 'animate', args: [[name],
+ {mode: 'immediate', frame: {duration: 0, redraw: true}, transition: {duration: 0}}]};
+ })}] : []
+ });
+ var pause = el._transitionData ? Plotly.animate(el, [null],
+ {mode: 'immediate', frame: {duration: 0}, transition: {duration: 0}}) : Promise.resolve();
+ pause.then(function() {
+ var previous = el._transitionData && el._transitionData._frames;
+ return previous && previous.length ? Plotly.deleteFrames(el) : null;
+ }).then(function() { return Plotly.react(el, cell.data, layout); })
+ .then(function() { return config.animated ? Plotly.addFrames(el, frames) : null; })
+ .then(function() {
+ busy = false;
+ if (pending) { pending = false; draw(); }
+ }).catch(function(error) { busy = false; console.warn(error); });
+ }
+ el.on('plotly_buttonclicked', function(event) {
+ if (!event.menu || event.menu.name !== 'tam-residual-type') return;
+ state.view = event.button.args[0];
+ draw();
+ });
+ el.on('plotly_animatingframe', function(event) {
+ var index = config.frames.indexOf(event.name);
+ if (index >= 0) {
+ state.frame = index;
+ if (busy) pending = true;
+ }
+ });
+ draw();
+}
diff --git a/inst/rmd/model_structure.Rmd b/inst/rmd/model_structure.Rmd
new file mode 100644
index 00000000..b46981ab
--- /dev/null
+++ b/inst/rmd/model_structure.Rmd
@@ -0,0 +1,8 @@
+Model structure
+================================================================================
+
+### Model structure
+
+```{r}
+plot_model(model_list = fits)
+```
diff --git a/inst/rmd/vis_tam.Rmd b/inst/rmd/vis_tam.Rmd
index 22fee565..80cf8c69 100644
--- a/inst/rmd/vis_tam.Rmd
+++ b/inst/rmd/vis_tam.Rmd
@@ -124,11 +124,7 @@ plot_mortality <- function(data, mean_line = list(dash = "dot", width = .5), ...
make_resid_plots <- function(d) {
fitted_models <- names(fits)[vapply(fits, inherits, logical(1), "tam_fit")]
d <- d[d$model %in% fitted_models, , drop = FALSE]
- if ("osa_res" %in% names(d) && !any(is.finite(d$osa_res))) {
- d$osa_res <- NULL
- }
- residual <- if ("osa_res" %in% names(d)) "osa_res" else "std_res"
- d <- droplevels(plot_rows(d, residual))
+ d <- droplevels(d)
if (!nrow(d)) {
return(setNames(rep(list(plotly::plot_ly()), 5),
c("bubbles", "year", "age", "cohort", "expected")))
@@ -178,6 +174,9 @@ table_link <- function(x, table_name = NULL, download_text = "Download table") {
```{r tables, child = "tables.Rmd"}
```
+```{r model-structure, child = "model_structure.Rmd"}
+```
+
diff --git a/man/add_osa_residuals.Rd b/man/add_osa_residuals.Rd
new file mode 100644
index 00000000..459fb5ea
--- /dev/null
+++ b/man/add_osa_residuals.Rd
@@ -0,0 +1,56 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/osa.R
+\name{add_osa_residuals}
+\alias{add_osa_residuals}
+\title{Calculate and retain one-step-ahead residuals}
+\usage{
+add_osa_residuals(fit, method = "oneStepGaussianOffMode", ...)
+}
+\arguments{
+\item{fit}{A fitted \code{tam_fit}.}
+
+\item{method}{One algorithm supported by \code{\link[RTMB:OSA-residuals]{RTMB::oneStepPredict()}}:
+\code{"oneStepGaussianOffMode"} (default), \code{"oneStepGaussian"},
+\code{"oneStepGeneric"}, \code{"fullGaussian"}, or \code{"cdf"}.}
+
+\item{...}{Named computation controls for \code{\link[RTMB:OSA-residuals]{RTMB::oneStepPredict()}}, such as
+\code{trace}, \code{parallel}, \code{subset}, \code{conditional}, or \code{splineApprox}.
+\code{subset} and \code{conditional} index included log observations in likelihood
+order. Their unrequested rows retain missing residuals.}
+}
+\value{
+A \code{tam_fit} with \code{osa_res} and \code{osa_method} columns in \verb{$obs_pred},
+and calculation metadata in \verb{$osa} and the observation-table attributes.
+}
+\description{
+Calculate observation diagnostics at fitted parameter values and return a fit
+containing the results. Use the returned object in \code{\link[=vis_tam]{vis_tam()}} or \code{\link[=residuals]{residuals()}}.
+No model is refitted, and the supplied fit is not modified.
+}
+\details{
+The algorithms use different predictive approximations and can
+differ numerically. \code{oneStepGeneric} uses numerical integration and can be
+expensive. The existing off-mode Gaussian calculation remains the default.
+Log observations have a continuous Gaussian likelihood, so \code{discrete = FALSE}
+is supplied unless explicitly overridden. Method-dependent upstream defaults,
+including \code{reverse}, are retained. When SAM is also loaded, parallel OSA
+calculations use the existing sequential safeguard.
+
+Temporary RTMB objectives protect the fitted parameters, uncertainty and
+objective state. Only OSA columns and metadata are added to a copy of the
+existing observation tables; row order, confidence intervals, units and custom
+columns are preserved. Metadata in \code{fit$osa} records the method, requested and
+effective settings, and observation coverage. Failed requested calculations
+remain missing and are recorded separately from unrequested observations.
+A calculation error returns no partially updated fit.
+}
+\examples{
+\dontrun{
+fit_osa <- add_osa_residuals(fit, method = "oneStepGeneric", trace = FALSE)
+residuals(fit_osa, type = "osa")
+vis_tam(fit_osa)
+}
+}
+\seealso{
+\code{\link[=residuals.tam_fit]{residuals.tam_fit()}}, \code{\link[=tidy_obs_pred]{tidy_obs_pred()}}
+}
diff --git a/man/check_obs.Rd b/man/check_obs.Rd
index 758c77d8..b6b550ca 100644
--- a/man/check_obs.Rd
+++ b/man/check_obs.Rd
@@ -4,10 +4,12 @@
\alias{check_obs}
\title{Check the structure of a TAM obs list}
\usage{
-check_obs(obs)
+check_obs(obs, quiet = FALSE)
}
\arguments{
\item{obs}{A named list of tinyAM input tables.}
+
+\item{quiet}{Logical. Suppress the success summary? Defaults to \code{FALSE}.}
}
\value{
\code{TRUE} (invisibly) if validation passes; otherwise aborts.
diff --git a/man/check_tam_data.Rd b/man/check_tam_data.Rd
index 29c8336f..ec35a2b5 100644
--- a/man/check_tam_data.Rd
+++ b/man/check_tam_data.Rd
@@ -4,11 +4,13 @@
\alias{check_tam_data}
\title{Check tinyAM input data}
\usage{
-check_tam_data(data)
+check_tam_data(data, quiet = FALSE)
}
\arguments{
\item{data}{A named list expected to contain data.frames \code{catch}, \code{index},
\code{weight}, and \code{maturity} with the columns described above.}
+
+\item{quiet}{Logical. Suppress the success summary? Defaults to \code{FALSE}.}
}
\value{
\code{TRUE} (invisibly) if validation passes; otherwise aborts.
@@ -47,11 +49,13 @@ requires nonempty labels on every row. Units never change input values.
See \code{\link[=as_tam_data]{as_tam_data()}} for units of derived quantities.
}
-On failure the function aborts with a \code{cli} error. On success, it returns
-\code{TRUE} (invisibly).
+On failure the function aborts with a \code{cli} error. On success, it prints a
+concise summary and returns \code{TRUE} invisibly. Use \code{quiet = TRUE} to suppress
+the success summary; validation errors are always reported.
}
\examples{
data(cod_obs)
-check_obs(cod_obs) # returns TRUE (invisibly) if valid
+check_tam_data(cod_obs)
+check_tam_data(cod_obs, quiet = TRUE)
}
diff --git a/man/fit_tam.Rd b/man/fit_tam.Rd
index 41d6a70f..888778b9 100644
--- a/man/fit_tam.Rd
+++ b/man/fit_tam.Rd
@@ -264,9 +264,9 @@ See \link{tinyAM-model} for the catch solver and feasibility safeguards.
\item{interval}{Level in \verb{(0, 1)} to use to generate confidence intervals,
where applicable; default \verb{0.95.}}
-\item{add_osa_res}{Logical; add one-step-ahead residuals? Hard wired to
-apply the \code{"oneStepGaussianOffMode"} method. See \code{\link[RTMB:OSA-residuals]{RTMB::oneStepPredict()}}
-for details.}
+\item{add_osa_res}{Logical; add one-step-ahead residuals using
+\code{"oneStepGaussianOffMode"}? To select another method after fitting and
+retain its results, use \code{\link[=add_osa_residuals]{add_osa_residuals()}}.}
\item{silent}{Logical; if \code{TRUE}, disables tracing information.}
diff --git a/man/get_population.Rd b/man/get_population.Rd
new file mode 100644
index 00000000..d83f68ea
--- /dev/null
+++ b/man/get_population.Rd
@@ -0,0 +1,73 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/getters.R
+\name{get_population}
+\alias{get_population}
+\alias{get_ssb}
+\alias{get_recruitment}
+\alias{get_fbar}
+\alias{get_mbar}
+\alias{get_f}
+\alias{get_m}
+\title{Extract common population estimates}
+\usage{
+get_ssb(..., model_list = NULL, interval = 0.95)
+
+get_recruitment(..., model_list = NULL, interval = 0.95)
+
+get_fbar(..., model_list = NULL, interval = 0.95)
+
+get_mbar(..., model_list = NULL, interval = 0.95)
+
+get_f(..., model_list = NULL, interval = 0.95)
+
+get_m(..., model_list = NULL, interval = 0.95)
+}
+\arguments{
+\item{...}{One or more fitted \code{tam_fit} models or \code{tam_ref} reporting
+references, or one named list of these objects.}
+
+\item{model_list}{A named list of models, instead of \code{...}.}
+
+\item{interval}{Confidence level in \verb{(0, 1)}, default \code{0.95}. Fitted summaries
+are rebuilt when their cached interval differs. Reference uncertainty is
+retained as supplied, because it cannot be recalculated from a reporting
+reference.}
+}
+\value{
+A data frame retaining available keys, \code{est}, confidence limits,
+\code{se}, \code{se_scale}, \code{unit} and \code{is_proj}. Available metadata, such as SSB
+timing, is retained. Annual summaries have \code{year}; surfaces also have \code{age}.
+}
+\description{
+Return a ready-to-plot table of spawning biomass, recruitment or mortality,
+without navigating the nested tidy output. These functions select the same
+estimates, uncertainty and units used in \code{\link[=tidy_tam]{tidy_tam()}} and the dashboard.
+}
+\details{
+\code{get_ssb()} and \code{get_recruitment()} return annual estimates. \code{get_fbar()}
+and \code{get_mbar()} return the model's abundance-weighted means over the ages
+chosen by \code{F_settings$mean_ages} and \code{M_settings$mean_ages}, respectively.
+\code{get_f()} and \code{get_m()} return year-by-age surfaces; F is total fishing
+mortality, summed over fleets. They do not return process deviations.
+
+A single fit returns its table without a model column. Multiple fits and
+named-list inputs add \code{model}, using the supplied names. Missing reference
+quantities are omitted rather than reconstructed. Explicit reference
+\code{comparison_scales} follow \code{\link[=tidy_tam]{tidy_tam()}}'s conversion rules; otherwise known
+incompatible units must be converted before combining models.
+
+These are lightweight selectors. They never refit models or calculate
+observation OSA residuals or joint-sample process diagnostics. Supplied
+objects are not modified.
+}
+\examples{
+\dontrun{
+get_ssb(fit)
+get_recruitment(fit, interval = 0.8)
+get_fbar(Base = fit, Alternative = other_fit)
+get_f(model_list = list(Base = fit, Alternative = other_fit))
+}
+}
+\seealso{
+\code{\link[=tidy_pop]{tidy_pop()}}, \code{\link[=tidy_tam]{tidy_tam()}}, \code{\link[=plot_trend]{plot_trend()}}, \code{\link[=plot_matrix]{plot_matrix()}}
+}
diff --git a/man/plot_model.Rd b/man/plot_model.Rd
new file mode 100644
index 00000000..b745775f
--- /dev/null
+++ b/man/plot_model.Rd
@@ -0,0 +1,51 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/model-diagram.R
+\name{plot_model}
+\alias{plot_model}
+\title{View the structure of fitted assessment models}
+\usage{
+plot_model(..., model_list = NULL)
+}
+\arguments{
+\item{...}{One or more \code{tam_fit} models or \code{tam_ref} reporting references,
+or one named list of these objects.}
+
+\item{model_list}{A named list of models, instead of \code{...}.}
+}
+\value{
+A browsable HTML container with a responsive SVG diagram and model
+settings. Display it directly, or save it with \code{\link[htmltools:save_html]{htmltools::save_html()}}.
+}
+\description{
+See how recruitment, survival, mortality and observations are connected,
+and which components have annual process variation. The diagram and its
+accompanying cards use the completed settings stored in each fitted model.
+}
+\details{
+Multiple models have a selector, initially showing the first tinyAM fit.
+A reporting reference displays an unavailable message because it does not
+contain a fitted tinyAM specification. The same viewer appears on the
+dashboard's Model structure page.
+
+An annual process set to \code{"off"} is labeled deterministic: it removes
+process deviations, not the underlying recruitment, survival or mortality.
+Initial states are shown separately from annual processes. Stock-recruit
+relationships add a parent-SSB feedback arrow. Formula cards retain the
+actual formulas, fleets, surveys, state blocks and projection controls.
+Short formula labels appear in the diagram; dashed input arrows identify
+the catch, index or population table supplying active predictors.
+
+This is a model overview, not a diagnostic of identifiability or numerical
+convergence. It reads stored settings without evaluating the fitted
+objective, sampling uncertainty, calculating residuals or changing the fit.
+No additional diagramming package is required.
+}
+\examples{
+\dontrun{
+plot_model(fit)
+plot_model(Base = fit, Alternative = other_fit)
+}
+}
+\seealso{
+\code{\link[=vis_tam]{vis_tam()}}, \code{\link[=fit_tam]{fit_tam()}}, \link{tinyAM-model}
+}
diff --git a/man/plot_tam.Rd b/man/plot_tam.Rd
index 9ab8e2d6..4da0b5e6 100644
--- a/man/plot_tam.Rd
+++ b/man/plot_tam.Rd
@@ -37,9 +37,23 @@ plot_obs_pred(
...
)
-plot_bubbles(data, res = NULL, y = ~age, ylab = "Age", title = NULL, ...)
+plot_bubbles(
+ data,
+ res = NULL,
+ y = ~age,
+ ylab = "Age",
+ title = NULL,
+ residual_type = c("auto", "standardized", "osa"),
+ ...
+)
-plot_resid(data, x = ~year, xlab = "Year", ...)
+plot_resid(
+ data,
+ x = ~year,
+ xlab = "Year",
+ residual_type = c("auto", "standardized", "osa"),
+ ...
+)
plot_par(data, ...)
}
@@ -70,6 +84,13 @@ applicable, Show CI / Hide CI.}
\item{res}{Optional one-column formula for signed bubble values, e.g.
\code{~est} for process deviations. Omitted values use OSA or standardized
observation residuals. Use \code{y = ~age_block} for shared M state blocks.}
+
+\item{residual_type}{Initial observation residual view: \code{"auto"} (default),
+\code{"standardized"} or \code{"osa"}. Auto uses available cached OSA residuals
+separately for each model in the supplied panel, otherwise standardized
+residuals. Partial OSA gaps remain missing. A dropdown switches views;
+unavailable OSA values are identified without computing them. A supplied
+\code{res} formula in \code{plot_bubbles()} bypasses these observation controls.}
}
\value{
A Plotly object.
diff --git a/man/residuals.tam_fit.Rd b/man/residuals.tam_fit.Rd
index 02ef92ce..e9064fdf 100644
--- a/man/residuals.tam_fit.Rd
+++ b/man/residuals.tam_fit.Rd
@@ -4,7 +4,12 @@
\alias{residuals.tam_fit}
\title{Extract observation residuals}
\usage{
-\method{residuals}{tam_fit}(object, type = c("standardized", "log", "response", "osa"), ...)
+\method{residuals}{tam_fit}(
+ object,
+ type = c("standardized", "log", "response", "osa"),
+ method = "oneStepGaussianOffMode",
+ ...
+)
}
\arguments{
\item{object}{A fitted \code{tam_fit}.}
@@ -12,17 +17,20 @@
\item{type}{\code{"standardized"} (default), \code{"log"} (log observation minus log
prediction), \code{"response"} (observation minus prediction), or \code{"osa"}.}
-\item{...}{Computation controls passed to \code{\link[=tidy_obs_pred]{tidy_obs_pred()}} and
-\code{\link[RTMB:OSA-residuals]{RTMB::oneStepPredict()}} when uncached OSA residuals are requested, such as
+\item{method}{OSA algorithm, as in \code{\link[=add_osa_residuals]{add_osa_residuals()}}. Only used with
+\code{type = "osa"}. Explicitly supplying it requests a fresh calculation.}
+
+\item{...}{Computation controls passed to \code{\link[RTMB:OSA-residuals]{RTMB::oneStepPredict()}}, such as
\code{trace = FALSE} or \code{parallel = FALSE}.
\code{subset} and \code{conditional} index included observations in likelihood order;
unrequested rows retain missing residuals.
-Controls are not used for cached OSA values; other residual types reject them.}
+Explicit controls request a fresh calculation; other residual types reject them.}
}
\value{
A data frame with year, age, type (catch/index), survey, observations,
-predictions, residuals, and available unit/fleet columns. The \code{residual_type}
-attribute names the requested scale.
+predictions, residuals, and available unit/fleet columns. OSA output also
+identifies \code{osa_method}. The \code{residual_type} attribute names the requested
+scale, and \code{osa} retains available method/settings/coverage metadata.
}
\description{
Return keyed catch and survey rows for plotting or further checks. This is a
@@ -30,20 +38,21 @@ convenience extractor: standardized residuals are the existing fitted-state
residuals, affected by shrinkage, not independent prediction errors.
}
\details{
-\code{type = "osa"} reuses cached OSA columns, or explicitly requests
-their calculation at fitted parameter values when absent. It can be expensive,
+\code{type = "osa"} reuses cached OSA columns when no computation controls
+or method are supplied. Otherwise both catch and index residuals are calculated
+together, using one algorithm at fitted parameter values. It can be expensive,
but does not rerun \code{\link[=fit_tam]{fit_tam()}} or re-estimate fixed parameters. Temporary
objectives protect the fitted object's state. Other residual types remain
simple extractors. Excluded observations, including zeros, retain rows with
missing residuals. The supplied fit is not modified or implicitly cached.
-To retain OSA results for later extraction and dashboard use, assign
-\code{fit$obs_pred <- tidy_obs_pred(fit, add_osa_res = TRUE, trace = FALSE)}.
+To retain results for later extraction and dashboard use, assign
+\code{fit_osa <- add_osa_residuals(fit, method = "oneStepGeneric", trace = FALSE)}.
}
\examples{
\dontrun{
residuals(fit)
residuals(fit, type = "log")
-residuals(fit, type = "osa", trace = FALSE)
+residuals(fit, type = "osa", method = "oneStepGeneric", trace = FALSE)
}
}
\seealso{
diff --git a/man/tidy_obs_pred.Rd b/man/tidy_obs_pred.Rd
index 23e446ee..b3d5b3fb 100644
--- a/man/tidy_obs_pred.Rd
+++ b/man/tidy_obs_pred.Rd
@@ -4,17 +4,24 @@
\alias{tidy_obs_pred}
\title{Tidy observed, predicted, and residual diagnostics}
\usage{
-tidy_obs_pred(fit, add_osa_res = FALSE, interval = 0.95, ...)
+tidy_obs_pred(
+ fit,
+ add_osa_res = FALSE,
+ interval = 0.95,
+ method = "oneStepGaussianOffMode",
+ ...
+)
}
\arguments{
\item{fit}{A fitted TAM object as returned by \code{\link[=fit_tam]{fit_tam()}}.}
-\item{add_osa_res}{Logical; add one-step-ahead residuals? Hard wired to
-apply the \code{"oneStepGaussianOffMode"} method.
-See \code{\link[RTMB:OSA-residuals]{RTMB::oneStepPredict()}} for details.}
+\item{add_osa_res}{Logical; calculate one-step-ahead residuals?}
\item{interval}{Confidence level for combined catchability intervals.}
+\item{method}{OSA algorithm; defaults to \code{"oneStepGaussianOffMode"}.
+See \code{\link[=add_osa_residuals]{add_osa_residuals()}} for supported choices and computation details.}
+
\item{...}{Arguments to pass to \code{\link[RTMB:OSA-residuals]{RTMB::oneStepPredict()}}.}
}
\value{
@@ -44,6 +51,10 @@ uncertainty in fitted states. Optional one-step-ahead residuals use a
different predictive calculation. See \link{tinyAM-model} for the observation equations.
OSA calculation rebuilds a temporary objective at stored fitted parameters;
it does not refit fixed parameters or modify the fitted objective.
+With \code{add_osa_res = TRUE}, it computes the requested method afresh and records
+its settings and observation coverage. \code{osa_method} identifies the algorithm
+beside \code{osa_res}. To retain diagnostics while preserving existing observation
+tables and confidence intervals, use \code{\link[=add_osa_residuals]{add_osa_residuals()}}.
When stockassessment is loaded, one-step residuals run sequentially: the
upstream TMB parallel helper cannot select between SAM and RTMB libraries.
}
diff --git a/man/vis_tam.Rd b/man/vis_tam.Rd
index cfc5c383..9eb4f72b 100644
--- a/man/vis_tam.Rd
+++ b/man/vis_tam.Rd
@@ -83,6 +83,9 @@ Supply models via \code{...} or \code{model_list}, but not both. The models must
form a uniquely named list so they can be labeled in the dashboard.
Additional arguments for \code{\link[rmarkdown:render]{rmarkdown::render()}} can be passed through
\code{render_args}, which must itself be a (named) list.
+Dashboard R packages are installed with tinyAM. Rendering also needs
+Pandoc, included with RStudio or available from
+\url{https://pandoc.org/installing.html}.
Precomputed assessment references may be mixed with tinyAM fits; they are
not refittable TAM objects. Supply \code{background} to add a page
describing assessment assumptions and translation choices.
@@ -93,6 +96,11 @@ only on the aggregate page; they are not allocated among fleets.
Without fleets, Aggregate fishery retains all single-fishery
diagnostics. With fleets, it shows total yield; residuals are
assessed separately within each fleet.
+Observation residuals default to Auto: each model uses retained OSA
+residuals when available in that panel, otherwise standardized
+residuals. The dropdown also selects either type explicitly. Partial
+OSA gaps remain missing, and labels identify the method. Retain OSA
+calculations with \code{\link[=add_osa_residuals]{add_osa_residuals()}}; rendering never computes them.
The Parameters menu separates Fixed and Random pages. Fixed plots
group similar quantities, with initial abundance shown separately.
Stock-recruit models add a relationship page with parent-year pairs,
diff --git a/tests/testthat/residual-selection-controls.js b/tests/testthat/residual-selection-controls.js
new file mode 100644
index 00000000..300a2173
--- /dev/null
+++ b/tests/testthat/residual-selection-controls.js
@@ -0,0 +1,75 @@
+const fs = require('fs'), assert = require('assert');
+const draw = eval('(' + fs.readFileSync(process.argv[2], 'utf8') + ')');
+const source = JSON.parse(fs.readFileSync(process.argv[3], 'utf8'));
+let banks = 0;
+global.Plotly = {
+ react: async (el, data, layout) => { el.data = data; el.layout = layout; },
+ addFrames: async (el, frames) => { banks++; el._transitionData = {_frames: frames}; },
+ deleteFrames: async el => { el._transitionData = {_frames: []}; },
+ animate: async () => {}
+};
+const flush = () => new Promise(resolve => setImmediate(resolve));
+function widget() {
+ const handlers = {};
+ return {handlers, on: (name, fun) => { handlers[name] = fun; }};
+}
+function chooseType(el, index) {
+ el.handlers.plotly_buttonclicked({menu: {name: 'tam-residual-type'}, button: {args: [index]}});
+}
+function chooseFrame(el, name) {
+ const frame = el._transitionData._frames.find(frame => frame.name === name);
+ el.data = frame.data;
+ Object.assign(el.layout, frame.layout);
+ el.handlers.plotly_animatingframe({name});
+}
+const points = el => el.data.reduce((n, trace) => n + trace.x.length, 0);
+(async () => {
+ const a = widget(), b = widget();
+ draw(a, {layout: {}}, JSON.parse(JSON.stringify(source)));
+ draw(b, {layout: {}}, JSON.parse(JSON.stringify(source)));
+ await flush();
+ assert.equal(points(a), 2);
+ assert.equal(points(b), 2);
+ assert.deepEqual(a.layout.updatemenus[0].buttons.map(button => button.label), ['Auto', 'Standardized', 'OSA']);
+ assert.equal(a.layout.updatemenus[0].x, 0);
+ assert.equal(a.layout.updatemenus[0].xanchor, 'left');
+ assert.equal(a.layout.updatemenus[0].yanchor, 'bottom');
+ assert.equal(a.layout.updatemenus[0].pad.b, 28);
+ assert.equal(a.layout.margin.t, 100); // Header rows stay clear of the upper-right modebar.
+ chooseFrame(a, 'B');
+ assert.equal(points(a), 6); // Auto keeps a model without cached OSA.
+ const before = banks;
+ chooseFrame(a, 'A');
+ chooseFrame(a, 'B');
+ assert.equal(banks, before); // Native frame changes keep their frame bank.
+ chooseType(a, 2); await flush();
+ assert.equal(a.layout.sliders[0].active, 1);
+ assert.equal(points(a), 0);
+ assert(a.layout.annotations[0].text.includes('OSA residuals unavailable'));
+ assert.equal(points(b), 2); // Widgets have independent selections.
+ chooseType(a, 1); await flush();
+ assert.equal(points(a), 6);
+ assert.equal(a.layout.sliders[0].active, 1);
+ chooseFrame(a, 'A');
+ assert.equal(points(a), 6);
+ chooseType(a, 0); await flush();
+ assert.equal(points(a), 2); // Partial OSA never fills its gaps from standard residuals.
+ assert(a.layout.annotations[0].text.includes('oneStepGeneric'));
+ for (let i = 0; i < 3; i++) chooseType(a, i);
+ await flush(); await flush();
+ assert.equal(a.layout.updatemenus[0].active, 2);
+ const singleConfig = JSON.parse(JSON.stringify(source));
+ singleConfig.frames = 'Only'; singleConfig.animated = false;
+ singleConfig.views.forEach(view => {
+ view.cells = [view.cells[0]];
+ view.cells[0].data.forEach(trace => {
+ trace.x = 2000; trace.y = 1; trace.text = 'one'; trace.customdata = 'custom';
+ trace.marker = trace.marker || {}; trace.marker.size = 1;
+ });
+ });
+ const single = widget(); draw(single, {layout: {}}, singleConfig); await flush();
+ assert.deepEqual(single.data[0].x, [2000]);
+ assert.deepEqual(single.data[0].customdata, ['custom']);
+ assert.deepEqual(single.layout.sliders, []);
+ console.log('Residual controls passed');
+})().catch(error => { console.error(error); process.exit(1); });
diff --git a/tests/testthat/test-fisheries-dashboard.R b/tests/testthat/test-fisheries-dashboard.R
index 6259ab1a..7d472670 100644
--- a/tests/testthat/test-fisheries-dashboard.R
+++ b/tests/testthat/test-fisheries-dashboard.R
@@ -67,7 +67,7 @@ test_that("Fisheries renders full fleet diagnostics and an unlabeled aggregate",
ref$pop$yield_pred_by_fleet$est <- ref$obs_pred$catch$pred
ref$fixed_par$est <- NA_real_
parameters <- default_fit$parameter_values
- reference_before <- serialize(ref,NULL)
+ reference_before <- ref
path <- tempfile(fileext=".html"); on.exit(unlink(path))
vis_tam(Unlabeled=default_fit,Fleet_reference=ref,output_file=path,
open_file=FALSE,render_args=list(quiet=TRUE))
@@ -92,5 +92,5 @@ test_that("Fisheries renders full fleet diagnostics and an unlabeled aggregate",
expect_match(aggregate_html,"Total yield",fixed=TRUE)
expect_false(grepl("Residuals (age, year)",aggregate_html,fixed=TRUE))
expect_identical(default_fit$parameter_values,parameters)
- expect_identical(serialize(ref,NULL),reference_before)
+ expect_identical(ref,reference_before)
})
diff --git a/tests/testthat/test-getters.R b/tests/testthat/test-getters.R
new file mode 100644
index 00000000..9ea9b155
--- /dev/null
+++ b/tests/testthat/test-getters.R
@@ -0,0 +1,79 @@
+test_that("quantity getters select the maintained annual and age-specific tables", {
+ getters <- list(ssb = get_ssb, recruitment = get_recruitment,
+ F_bar = get_fbar, M_bar = get_mbar, F = get_f, M = get_m)
+ original <- default_fit
+ expected <- .population_tables(list(base = default_fit), .95)[[1L]]
+ for (name in names(getters)) {
+ value <- getters[[name]](default_fit)
+ reference <- expected[[name]]
+ rownames(reference) <- NULL
+ expect_equal(value, reference)
+ expect_false("model" %in% names(value))
+ expect_true(all(c("year", "est", "is_proj") %in% names(value)))
+ }
+ expect_true("age" %in% names(get_f(default_fit)))
+ expect_false("age" %in% names(get_fbar(default_fit)))
+ expect_identical(default_fit, original)
+})
+
+test_that("quantity getters retain labels, intervals and projections without diagnostics", {
+ testthat::local_mocked_bindings(
+ tidy_obs_pred = function(...) stop("Unexpected observation diagnostics"),
+ tidy_par = function(...) stop("Unexpected parameter tidying"),
+ tidy_recruitment = function(...) stop("Unexpected recruitment diagnostics"),
+ .tidy_process_residuals = function(...) stop("Unexpected process diagnostics"))
+ original <- N_dev
+ value <- get_ssb(Base = default_fit, Forecast = N_dev, interval = .8)
+ expect_identical(names(value)[1L], "model")
+ expect_identical(unique(value$model), c("Base", "Forecast"))
+ expect_true(any(value$is_proj))
+ expected <- tidy_pop(N_dev, interval = .8)$ssb
+ expect_equal(value$upr[value$model == "Forecast"], expected$upr)
+ expect_identical(get_ssb(list(Base = default_fit))$model,
+ rep("Base", nrow(default_fit$pop$ssb)))
+ expect_identical(get_ssb(model_list = list(`2020` = default_fit))$model,
+ rep(2020L, nrow(default_fit$pop$ssb)))
+ expect_identical(N_dev, original)
+})
+
+test_that("population selectors share explicit reference conversions and SE scales", {
+ fit <- default_fit
+ fit$pop$ssb$unit <- "tonnes"
+ reference <- structure(list(pop = list(ssb = fit$pop$ssb),
+ comparison_scales = c(ssb = 1e-3)), class = c("tam_ref", "list"))
+ reference$pop$ssb$unit <- "thousand tonnes"
+ original <- list(fit = fit, reference = reference)
+ value <- get_ssb(model_list = list(Base = fit, Accepted = reference))
+ tiny <- value[value$model == "Base", ]
+ expect_equal(tiny$est, fit$pop$ssb$est * 1e-3)
+ expect_equal(tiny$se, fit$pop$ssb$se)
+ expect_true(all(value$unit == "thousand tonnes"))
+ fit$pop$ssb$se_scale <- "reported"
+ expect_equal(get_ssb(fit, reference)$se[seq_len(nrow(fit$pop$ssb))],
+ fit$pop$ssb$se * 1e-3)
+ expect_identical(reference, original$reference)
+ bad <- reference
+ bad$comparison_scales <- NULL
+ expect_error(get_ssb(default_fit, bad), "incompatible units")
+ bad <- reference
+ bad$comparison_scales <- c(ssb = .01)
+ expect_error(get_ssb(default_fit, bad), "require")
+ missing <- structure(list(pop = list()), class = "tam_ref")
+ expect_equal(get_f(default_fit, missing)$est, get_f(default_fit)$est)
+ expect_error(get_f(missing), "unavailable")
+})
+
+test_that("quantity getters preserve configured means and missing uncertainty", {
+ fit <- N_dev
+ expect_equal(get_fbar(fit)$est, fit$pop$F_bar$est)
+ expect_equal(get_mbar(fit)$est, fit$pop$M_bar$est)
+ fit$pop <- NULL
+ fit$sdrep <- NULL
+ value <- get_ssb(fit)
+ expect_true(all(is.na(value$lwr)))
+ expect_true(all(is.na(value$upr)))
+ expect_equal(value$est, unname(fit$rep$ssb))
+ expect_error(get_ssb(fit, interval = 1), "strictly between")
+ expect_error(get_ssb(fit, model_list = list(Base = fit)), "either")
+ expect_error(get_ssb(model_list = list(fit)), "named")
+})
diff --git a/tests/testthat/test-model-diagram.R b/tests/testthat/test-model-diagram.R
new file mode 100644
index 00000000..26e8fe45
--- /dev/null
+++ b/tests/testthat/test-model-diagram.R
@@ -0,0 +1,101 @@
+test_that("model diagrams read completed settings without calculations or mutation", {
+ testthat::local_mocked_bindings(
+ fit_tam = function(...) stop("Unexpected refit"),
+ tidy_tam = function(...) stop("Unexpected diagnostics"),
+ nll_fun = function(...) stop("Unexpected objective evaluation"))
+ fit <- default_fit
+ original <- fit
+ set.seed(251)
+ seed <- .Random.seed
+ view <- plot_model(fit)
+ expect_s3_class(view, "shiny.tag")
+ html <- htmltools::renderTags(view)$html
+ expect_match(html, "Stochastic RW")
+ expect_match(html, "Deterministic: process off")
+ expect_match(html, "Older-age initialization: exp")
+ expect_match(html, "Recruitment boundary: free")
+ expect_match(html, "F is the sum of fleet F")
+ expect_match(html, "tinyAM-model.html", fixed = TRUE)
+ expect_match(html, "q formula:")
+ expect_false(grepl("Parent SSB:", html, fixed = TRUE))
+ expect_match(html, 'data-formula="survey"', fixed = TRUE)
+ expect_match(html, 'data-input-component="survey"', fixed = TRUE)
+ expect_match(html, 'index: q_block', fixed = TRUE)
+ expect_identical(fit, original)
+ expect_identical(.Random.seed, seed)
+})
+
+test_that("model diagrams show deterministic means and stock-recruit feedback", {
+ dat <- make_test_dat(years = 2000:2007, ages = 2:8,
+ N_settings = list(process = "off", init = "latent"),
+ rec_settings = list(process = "off", mu_form = ~ bh(ssb)),
+ F_settings = list(process = "off", mu_form = ~ factor(year) + factor(age)),
+ M_settings = list(process = "off", mu_supplied = ~ I(.3)))
+ fit <- structure(list(dat = dat), class = "tam_fit")
+ html <- htmltools::renderTags(plot_model(fit))$html
+ expect_match(html, "Mean: ~factor\\(year\\) \\+ factor\\(age\\)")
+ expect_match(html, "Parent SSB: bh, lag 2", fixed = TRUE)
+ expect_match(html, "Parent-SSB relationship: bh ; lag 2", fixed = TRUE)
+ expect_match(html, "Older-age initialization: latent", fixed = TRUE)
+ expect_match(html, "Deterministic (process off)", fixed = TRUE)
+ expect_false(grepl("Stochastic process:", html, fixed = TRUE))
+ expect_match(html, 'data-input-component="fishing"', fixed = TRUE)
+ expect_match(html, "catch: year, age", fixed = TRUE)
+ spec <- .model_structure_spec(dat)
+ expect_identical(unname(spec$inputs["recruitment"]), "")
+ dat$rec_settings$mu_form <- ~ factor(year)
+ expect_identical(unname(.model_structure_spec(dat)$inputs["recruitment"]), "")
+ dat$pop_data$temperature <- seq_len(nrow(dat$pop_data))
+ dat$N_settings$mu_form <- ~ temperature
+ dat$M_settings$mu_form <- ~ temperature
+ spec <- .model_structure_spec(dat)
+ expect_equal(unname(spec$inputs[c("population", "natural")]),
+ rep("pop_covariates: temperature", 2L))
+ expect_equal(unname(spec$formula_labels["natural"]), "Mean effects ~temperature")
+})
+
+test_that("model diagrams retain fleet, survey and shared block definitions", {
+ data <- cod_obs
+ data$catch <- rbind(transform(data$catch, fleet = "Trawl", obs = .6 * obs),
+ transform(data$catch, fleet = "Gillnet", obs = .4 * obs))
+ data$index <- rbind(transform(data$index, survey = "RV"),
+ transform(data$index, survey = "Sentinel"))
+ dat <- make_test_dat(data = data, years = 2000:2007, ages = 2:8,
+ F_settings = list(process = "off", mu_form = ~ 0 +
+ for_fleet("Trawl", ~ factor(year) + logistic(age)) +
+ for_fleet("Gillnet", ~ factor(year) + double_logistic(age))),
+ M_settings = list(process = "ar1", age_breaks = c(2, 5, 9)))
+ fit <- structure(list(dat = dat), class = "tam_fit")
+ html <- htmltools::renderTags(plot_model(fit))$html
+ expect_match(html, "Fleets: Gillnet, Trawl", fixed = TRUE)
+ expect_match(html, "Surveys: RV, Sentinel", fixed = TRUE)
+ expect_match(html, "for_fleet", fixed = TRUE)
+ expect_match(html, "double_logistic", fixed = TRUE)
+ expect_match(html, "Age blocks: 2-4, 5-8", fixed = TRUE)
+ expect_match(html, "Shared age blocks share actual M", fixed = TRUE)
+})
+
+test_that("model selectors handle references and independent escaped containers", {
+ reference <- structure(list(pop = list()), class = "tam_ref")
+ fits <- list(Accepted = reference, `Fit \"quoted\"` = default_fit)
+ html <- htmltools::renderTags(plot_model(model_list = fits))$html
+ expect_match(html, '