diff --git a/NAMESPACE b/NAMESPACE index f7fd3c70..cdeda6cd 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -17,11 +17,13 @@ export(plot_catch_comp) export(plot_discard) export(plot_fishing_mortality) export(plot_index) +export(plot_joint_residuals) export(plot_landings) export(plot_natural_mortality) export(plot_obsvpred) export(plot_recruitment) export(plot_recruitment_deviations) +export(plot_runs_test) export(plot_selectivity) export(plot_spawning_biomass) export(plot_spr) diff --git a/R/conout_fims.R b/R/conout_fims.R new file mode 100644 index 00000000..490000b1 --- /dev/null +++ b/R/conout_fims.R @@ -0,0 +1,143 @@ +#' Convert FIMS output to the standard output format +#' +#' Reshape the estimates from a FIMS model (`FIMS::get_estimates()`) into the +#' standardized long format returned by [convert_output()]. +#' +#' @param file A tibble returned by `FIMS::get_estimates()`, a `FIMSFit` +#' object, or the path to an .rds file containing either. +#' @param out_new The empty standard output data frame from +#' [convert_output()]. +#' @param start_year Number. First year of the model. FIMS reports years as +#' indices (`year_i`, 1 = first model year), which are converted to years as +#' `year_i + start_year - 1`. NULL keeps the indices. +#' +#' @returns A data frame with the columns of `out_new`. +#' +#' @details Each data point (index, landings, and age and length compositions) +#' is one row in FIMS, with the observed value in `observed` and the expected +#' value in `estimated`. These are split into two rows, labeled +#' "*_observed" and "*_predicted" as for other models (e.g., "index_observed" +#' and "index_predicted"). Observed rows hold the observation error in +#' `uncertainty` and the negative log-likelihood of the observation in +#' `likelihood`. FIMS reports the observation error as `log_sd`, the log of +#' its standard deviation. For lognormal data this standard deviation +#' (\eqn{\sigma}) is on the log scale, so it is converted to the coefficient +#' of variation of the observation, \eqn{\sqrt{\exp(\sigma^2) - 1}}, and +#' labeled "cv"; for other data `exp(log_sd)` is labeled "se". Observations +#' coded as missing (-999) are not observations and are not included. +#' Compositions are kept as counts, as reported by FIMS. +#' +#' All other estimates keep one row each, with the value in `estimate` and +#' its standard error in `uncertainty`. FIMS codes estimates that could not +#' be calculated as -999; these are NA. Rows in years after the last year of +#' data (e.g., spawning biomass at the start of the year after the model +#' ends) have `era` "fore"; other rows with a year have `era` "time". +#' +#' @noRd +conout_fims <- function(file, out_new, start_year = NULL) { + dat <- read_fims_estimates(file) + names(dat) <- tolower(names(dat)) + if (!all(c("label", "estimated") %in% names(dat))) { + cli::cli_abort(c( + "{.arg file} is not FIMS output.", + "i" = "Supply the output of {.code FIMS::get_estimates()} or a {.cls FIMSFit} object." + )) + } + col <- function(name) { + if (name %in% names(dat)) dat[[name]] else rep(NA, nrow(dat)) + } + year_i <- suppressWarnings(as.numeric(col("year_i"))) + if (is.null(start_year)) { + cli::cli_alert_info( + "FIMS reports years as indices (1 = first model year); use {.arg start_year} to convert them to years." + ) + year <- year_i + } else { + year <- year_i + start_year - 1 + } + observed <- suppressWarnings(as.numeric(col("observed"))) + observed[observed %in% -999] <- NA + # The last year of data; later years are after the model ends + data_years <- year_i[!is.na(observed)] + last_year <- if (length(data_years) > 0) max(data_years) else max(year_i, na.rm = TRUE) + era <- dplyr::case_when( + is.na(year_i) ~ NA_character_, + year_i > last_year ~ "fore", + TRUE ~ "time" + ) + uncertainty <- dplyr::na_if(suppressWarnings(as.numeric(col("uncertainty"))), -999) + base <- data.frame( + module_name = as.character(col("module_name")), + label = as.character(dat$label), + era = era, + year = year, + initial = suppressWarnings(as.numeric(col("input"))), + estimate = dplyr::na_if(suppressWarnings(as.numeric(dat$estimated)), -999), + uncertainty = uncertainty, + uncertainty_label = ifelse(is.na(uncertainty), NA_character_, "se"), + likelihood = NA_real_, + fleet = as.character(col("fleet")), + age = suppressWarnings(as.numeric(col("age"))), + len_bins = suppressWarnings(as.numeric(col("length"))), + stringsAsFactors = FALSE + ) + + # Data points: split into observed and predicted rows + stream <- dplyr::case_when( + grepl("^index_?expected$", base$label, ignore.case = TRUE) ~ "index", + grepl("^(catch|landings)_?expected$", base$label, ignore.case = TRUE) ~ "landings", + grepl("^age_?comp(osition)?_?expected$", base$label, ignore.case = TRUE) ~ "age_composition", + grepl("^length_?comp(osition)?_?expected$", base$label, ignore.case = TRUE) ~ "length_composition", + TRUE ~ NA_character_ + ) + # Fleets with data for a stream; expected values for other fleets (e.g., + # the expected index of a fishery without an index) keep their FIMS label + has_data <- stats::ave(!is.na(observed), stream, base$fleet, FUN = any) + is_data <- !is.na(stream) & has_data %in% TRUE + obs_rows <- base[is_data & !is.na(observed), , drop = FALSE] + if (nrow(obs_rows) > 0) { + is_obs <- is_data & !is.na(observed) + sigma <- exp(suppressWarnings(as.numeric(col("log_sd")))[is_obs]) + # A lognormal sigma is on the log scale; report the CV of the observation + lognormal <- grepl("^(d?lnorm|log_?normal)$", col("distribution")[is_obs], ignore.case = TRUE) + obs_rows$label <- paste0(stream[is_obs], "_observed") + obs_rows$estimate <- observed[is_obs] + obs_rows$uncertainty <- ifelse(lognormal, sqrt(exp(sigma^2) - 1), sigma) + obs_rows$uncertainty_label <- dplyr::case_when( + is.na(sigma) ~ NA_character_, + lognormal ~ "cv", + TRUE ~ "se" + ) + obs_rows$initial <- NA_real_ + # Stock Synthesis reports the negative log-likelihood of each observation + obs_rows$likelihood <- -suppressWarnings(as.numeric(col("likelihood")))[is_obs] + } + pred_rows <- base[is_data, , drop = FALSE] + pred_rows$label <- paste0(stream[is_data], "_predicted") + out <- rbind(base[!is_data, , drop = FALSE], obs_rows, pred_rows) + + missing <- setdiff(names(out_new), names(out)) + out[missing] <- NA + out[, names(out_new)] +} + +# The estimates of a FIMS model from get_estimates() output, a FIMSFit +# object, or an .rds file with either +read_fims_estimates <- function(file) { + if (is.character(file) && length(file) == 1) { + file <- readRDS(file) + } + if (inherits(file, "FIMSFit")) { + get_estimates <- tryCatch( + getExportedValue("FIMS", "get_estimates"), + error = function(e) { + cli::cli_abort("The {.pkg FIMS} package is needed to read a {.cls FIMSFit} object.") + } + ) + file <- get_estimates(file) + } + if (!is.data.frame(file)) { + cli::cli_abort("{.arg file} must be FIMS output from {.code FIMS::get_estimates()} or a {.cls FIMSFit} object.") + } + as.data.frame(file) +} diff --git a/R/convert_output.R b/R/convert_output.R index 62bda5b8..d3b2ac99 100644 --- a/R/convert_output.R +++ b/R/convert_output.R @@ -3,7 +3,8 @@ #' Format stock assessment output file to a standardized tibble. #' #' @param file Path. Assessment model output file path. Currently accepts SS3 -#' `r4ss::ss_output()` or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object. +#' `r4ss::ss_output()` or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output +#' (the tibble returned by `FIMS::get_estimates()`, a `FIMSFit` object, or an .rds file with either). #' @param model String. Assessment model used in evaluation #' #' Default: NULL @@ -20,6 +21,11 @@ #' in the file path provided. #' #' Default: NULL +#' @param start_year Number. First year of the model, for FIMS output only. +#' FIMS reports years as indices (1 = first model year), which are converted +#' to years when `start_year` is given. +#' +#' Default: NULL (years are left as indices) #' #' @author Samantha Bredeck #' @@ -30,6 +36,14 @@ #' version of a model output file. Converted data frame is always returned. #' It will also be saved if save_dir is not NULL. #' +#' For FIMS, each data point (index, landings, and age and length +#' compositions) is split into an observed row (e.g., "index_observed") and a +#' predicted row (e.g., "index_predicted"), as for the other models. Observed +#' rows hold the observation error (a log-scale standard error for lognormal +#' data) in `uncertainty` and the negative log-likelihood in `likelihood`. +#' Other estimates are relabeled with `fims_var_names.csv` (e.g., +#' "expected_recruitment" becomes "recruitment"). +#' #' #' @export #' @@ -46,7 +60,8 @@ convert_output <- function( file, model = NULL, fleet_names = NULL, - save_dir = NULL + save_dir = NULL, + start_year = NULL ) { if (is.character(file) && grepl("Report.sso", file)) { @@ -163,8 +178,12 @@ convert_output <- function( cli::cli_abort("Unknown file type. Please indicate model.") ) } - } else { + } else if (is.null(model)) { model <- switch(class(file)[1], + "FIMSFit" = { + cli::cli_alert_info("Processing FIMS output file...") + "fims" + }, "fims" = { cli::cli_alert_info("Processing FIMS output file...") "fims" @@ -210,17 +229,7 @@ convert_output <- function( cli::cli_abort("JABBA output not currently compatible.") #### FIMS #### } else if (model == "fims") { - if (is.character(file) && length(file) == 1) { - fims_output <- readRDS(file) - } else { - fims_output <- file - } - # TEMPORARY -- add in era to df - if ("era" %notin% colnames(fims_output)) { - fims_output$era <- NA - } - fims_output[setdiff(tolower(names(out_new)), tolower(names(fims_output)))] <- NA - out_new <- fims_output + out_new <- conout_fims(file = file, out_new = out_new, start_year = start_year) #### Rceattle #### } else if (model == "rceattle") { out_new <- conout_rceattle(file = file, out_new = out_new, fleet_names = fleet_names) diff --git a/R/plot_index.R b/R/plot_index.R index 4acb0a13..03552cbc 100644 --- a/R/plot_index.R +++ b/R/plot_index.R @@ -59,7 +59,9 @@ plot_index <- function( # Filter data prepared_data <- filter_data( dat, - label_name = "index", + # Only the observed and predicted index are plotted; other labels that + # contain "index" (e.g., catchability or index numbers at age) are not + label_name = "^index_(observed|predicted)$", era = NULL, geom = "line", # ifelse guarantees the code doesn't miss grouping when label has > 1 value @@ -90,6 +92,10 @@ plot_index <- function( if (!is.null(group)) { facet <- c(facet, group) } + # process_data() drops columns with a single value (e.g., an index from one + # fleet), so facet only by the columns that remain + facet <- intersect(facet, colnames(prepared_data)) + if (length(facet) == 0) facet <- NULL # transform prep data so group_var = label # prepared_data <- prepared_data |> # dplyr::mutate( @@ -117,7 +123,7 @@ plot_index <- function( ) ) # Overwrite facets from base plot_timeseries bc scales need to be free - if ("fleet" %in% colnames(prepared_data)) { + if (!is.null(facet) && "fleet" %in% colnames(prepared_data)) { facet <- paste("~", paste(facet, collapse = " + ")) facet_formula <- stats::reformulate(facet) plt <- plt + ggplot2::facet_wrap(facet_formula, scales = "free", labeller = function(labels) { diff --git a/R/plot_joint_residuals.R b/R/plot_joint_residuals.R new file mode 100644 index 00000000..aa1def07 --- /dev/null +++ b/R/plot_joint_residuals.R @@ -0,0 +1,329 @@ +#' Plot joint residuals of the fits to indices of abundance +#' +#' Plot the residuals of the fits to all indices over time with a loess +#' smoother and a fit statistic, as in `JABBA` (`jbplot_residuals()` and +#' `jbplot_stdresiduals()`) and `ss3diags` (`SSplotJABBAres()`). A trend in +#' the smoother shows a pattern in the residuals; a large statistic shows a +#' poor fit. +#' +#' @inheritParams plot_runs_test +#' @param scaled Logical. TRUE/FALSE; plot standardized residuals (log +#' residuals divided by their standard error) instead of log residuals. +#' +#' Default: FALSE +#' @param overlay Logical. TRUE/FALSE; plot all indices on one panel, with a +#' boxplot of the residuals in each year. FALSE plots each index on its own +#' panel. +#' +#' Default: TRUE +#' @param show_stat Logical. TRUE/FALSE; show the fit statistic in the top +#' right corner of each panel: the root mean square error (RMSE) of log +#' residuals, or the standard deviation of normalized residuals (SDNR) of +#' scaled residuals. When `overlay = TRUE`, the statistic of each index is +#' also shown in the legend. +#' +#' Default: TRUE +#' @param ncol Number. Number of columns of panels when `overlay = FALSE`. +#' +#' Default: NULL +#' +#' @returns A `ggplot2` object. +#' +#' @details Log residuals are \eqn{r = \log(o) - \log(p)} for observed +#' \eqn{o} and predicted \eqn{p}, as in [plot_runs_test()]; scaled +#' residuals are \eqn{z = r / s} for the log-scale standard error \eqn{s} of +#' the observation. The RMSE is \eqn{100 \sqrt{\overline{r^2}}} percent. The +#' SDNR is \eqn{\sqrt{\sum z^2 / (n - 1)}}, shown with its critical value +#' \eqn{\sqrt{\chi^2_{0.95, n - 1} / (n - 1)}}; an SDNR above the critical +#' value suggests the standard errors are too small (Francis 2011). +#' +#' The smoother is a loess regression (span 0.75) of the mean residual of +#' each year on year, for all indices together or for each index on its own +#' panel. It is not drawn for fewer than five years of residuals, or when +#' the years are too sparse for a stable fit; a message names the indices +#' without one. +#' +#' Observations are treated as in [plot_runs_test()]: a fleet observed in +#' more than one area, season, subseason, month, or sex becomes separate +#' series, and observations without a positive observed and predicted value +#' are reported in a message. Scaled residuals also need a standard error; +#' observations without one are reported in a message. +#' +#' @references Francis, R.I.C.C. 2011. Data weighting in statistical +#' fisheries stock assessment models. Canadian Journal of Fisheries and +#' Aquatic Sciences 68: 1124-1138. +#' +#' @seealso [plot_runs_test()], [plot_index()], [convert_output()] +#' @export +#' +#' @examples +#' plot_joint_residuals(stockplotr::example_data, interactive = FALSE) +#' plot_joint_residuals( +#' stockplotr::example_data, +#' scaled = TRUE, +#' overlay = FALSE, +#' interactive = FALSE +#' ) +plot_joint_residuals <- function( + dat, + fleet = NULL, + module = NULL, + scaled = FALSE, + overlay = TRUE, + show_stat = TRUE, + ncol = NULL, + interactive = TRUE, + make_rda = FALSE, + figures_dir = getwd() +) { + check_data(dat) + resid <- index_residuals(dat, fleet = fleet, module = module, interactive = interactive) + if (scaled) { + no_se <- is.na(resid$se) | resid$se <= 0 + if (any(no_se)) { + cli::cli_alert_info( + "{sum(no_se)} residual{?s} without a standard error {?is/are} not shown when {.arg scaled} is TRUE." + ) + resid <- resid[!no_se, , drop = FALSE] + if (nrow(resid) == 0) { + cli::cli_abort("No residuals have a standard error to scale by.") + } + } + resid$residual <- resid$residual / resid$se + } + + stat <- residual_stat(resid, scaled = scaled) + by_series <- stat[stat$series != "All indices", , drop = FALSE] + resid$index <- factor( + resid$series, + levels = by_series$series, + # The legend shows each index's statistic; panels show it in the corner + labels = if (show_stat && overlay) { + paste0(by_series$series, " (", by_series$value, ")") + } else { + by_series$series + } + ) + smooth <- residual_smooth(resid, by_series = !overlay) + + plt <- ggplot2::ggplot(resid, ggplot2::aes(x = .data$year, y = .data$residual)) + + ggplot2::geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") + if (overlay) { + n_series <- nrow(by_series) + dodge <- ggplot2::position_dodge(width = if (n_series > 1) 0.6 else 0) + # Boxplots only where a year has more than one residual + boxes <- resid[resid$year %in% resid$year[duplicated(resid$year)], , drop = FALSE] + plt <- plt + + ggplot2::geom_boxplot( + data = boxes, + ggplot2::aes(group = .data$year), + fill = "gray85", + color = "gray60", + outliers = FALSE, + width = 0.8 + ) + } else { + dodge <- "identity" + } + plt <- plt + + ggplot2::geom_linerange( + ggplot2::aes(ymin = 0, ymax = .data$residual, color = .data$index), + position = dodge + ) + + ggplot2::geom_point( + ggplot2::aes(fill = .data$index), + shape = 21, + size = 2, + position = dodge + ) + if (!is.null(smooth)) { + plt <- plt + + ggplot2::geom_line( + data = smooth, + ggplot2::aes(x = .data$year, y = .data$smooth), + linewidth = 1 + ) + } + if (show_stat) { + labels <- if (overlay) { + stat[stat$series == "All indices", , drop = FALSE] + } else { + by_series + } + labels$index <- if (overlay) NA else resid$index[match(labels$series, resid$series)] + plt <- plt + + ggplot2::geom_label( + data = labels, + ggplot2::aes(x = Inf, y = Inf, label = .data$text), + hjust = 1.05, vjust = 1.3, + fill = "white", alpha = 0.8, label.size = 0, + inherit.aes = FALSE + ) + } + legend_ncol <- if (max(nchar(levels(resid$index))) > 30) 1 else 2 + plt <- plt + + ggplot2::labs( + x = "Year", + y = if (scaled) "Standardized residual" else "Residual (log observed / predicted)", + color = "Index", + fill = "Index" + ) + + ggplot2::guides( + color = ggplot2::guide_legend(ncol = legend_ncol), + fill = ggplot2::guide_legend(ncol = legend_ncol) + ) + + theme_noaa() + + ggplot2::theme(legend.position = "bottom") + if (!overlay) { + plt <- plt + + ggplot2::facet_wrap(ggplot2::vars(.data$index), ncol = ncol) + + ggplot2::theme(legend.position = "none") + } + + if (make_rda) { + fit_export_rda( + plt, + topic_label = "joint_residuals", + tags = joint_residuals_tags(resid, stat, smooth, scaled = scaled, overlay = overlay), + kqs = joint_residuals_kqs(resid, stat, scaled = scaled), + dat = dat, + figures_dir = figures_dir + ) + } + plt +} + +# Fit statistic of the residuals of each series and of all residuals +# together ("All indices", the last row): the RMSE (percent) of log +# residuals, as in ss3diags::SSrmse(), or the SDNR of scaled residuals, with +# its critical value, as in JABBA::jbplot_stdresiduals() +residual_stat <- function(resid, scaled = FALSE) { + stat_fun <- function(r) { + n <- length(r) + if (scaled) { + sdnr <- if (n > 1) sqrt(sum(r^2) / (n - 1)) else NA_real_ + critical <- if (n > 1) sqrt(stats::qchisq(0.95, df = n - 1) / (n - 1)) else NA_real_ + tibble::tibble(stat = sdnr, critical = critical, n = n) + } else { + tibble::tibble(stat = 100 * sqrt(mean(r^2)), critical = NA_real_, n = n) + } + } + series <- unique(resid$series) + out <- dplyr::bind_rows( + lapply(series, function(s) stat_fun(resid$residual[resid$series == s])), + stat_fun(resid$residual) + ) + out <- tibble::tibble(series = c(series, "All indices"), out) + if (scaled) { + out$value <- ifelse( + is.na(out$stat), + "NA", + paste0(format(round(out$stat, 2), nsmall = 2, trim = TRUE)) + ) + out$text <- paste0( + "SDNR = ", out$value, + ifelse(is.na(out$critical), "", paste0(" (critical ", format(round(out$critical, 2), nsmall = 2, trim = TRUE), ")")) + ) + } else { + out$value <- paste0(format(round(out$stat, 1), nsmall = 1, trim = TRUE), "%") + out$text <- paste0("RMSE = ", out$value) + } + out +} + +# Loess smoother of the mean residual of each year, for all residuals or for +# each series; NULL (or no rows for a series) for fewer than five years +residual_smooth <- function(resid, by_series = FALSE) { + groups <- if (by_series) c("index", "series") else character(0) + yearly <- resid |> + dplyr::summarise(mean = mean(.data$residual), .by = dplyr::all_of(c(groups, "year"))) |> + dplyr::arrange(.data$year) + # NULL for fewer than five years, or when loess cannot fit them stably + # (it warns when the years are too sparse for the span) + fit_one <- function(d) { + if (nrow(d) < 5) { + return(NULL) + } + tryCatch( + { + fit <- stats::loess(mean ~ year, data = d, span = 0.75) + d$smooth <- stats::predict(fit, newdata = d) + d + }, + warning = function(w) NULL + ) + } + parts <- if (by_series) lapply(split(yearly, yearly$series), fit_one) else list(fit_one(yearly)) + short <- if (by_series) setdiff(unique(yearly$series), names(Filter(Negate(is.null), parts))) else NULL + out <- dplyr::bind_rows(parts) + if (nrow(out) == 0) { + cli::cli_alert_info("Too few years of residuals for a loess smoother; it is not drawn.") + return(NULL) + } + if (length(short) > 0) { + cli::cli_alert_info("Too few years of residuals for a loess smoother for {short}; none is drawn for {?it/them}.") + } + out +} + +# Tags for the caption and alternative text (joint.residuals.*) +joint_residuals_tags <- function(resid, stat, smooth, scaled, overlay) { + series <- stat[stat$series != "All indices", , drop = FALSE] + all <- stat[stat$series == "All indices", , drop = FALSE] + worst <- series[which.max(series$stat), , drop = FALSE] + stat_name <- if (scaled) { + "standard deviation of normalized residuals (SDNR)" + } else { + "root mean square error (RMSE)" + } + trend <- if (is.null(smooth)) { + "There are too few years to draw a smoother." + } else { + range <- range(smooth$smooth) + paste0( + "The smoother through the mean residual of each year ranges from ", + round(range[1], 2), " to ", round(range[2], 2), "." + ) + } + tags <- list( + type = if (scaled) "Standardized residuals" else "Residuals", + index.list = fit_list_text(series$series), + layout = paste0( + if (overlay) { + "together over time, with a boxplot of the residuals in each year (gray) and a loess smoother through the mean residual of each year (black line)" + } else { + "over time, each index on its own panel with a loess smoother (black line)" + }, + ". Residuals are the log of the observed index divided by the predicted index", + if (scaled) ", divided by the standard error of the observation" else "", + ", and the ", stat_name, " of each index is ", + if (overlay) "in the legend" else "in the corner of its panel", + if (scaled) "; an SDNR above its critical value suggests the standard errors are too small" else "" + ), + start.year = min(resid$year), + end.year = max(resid$year), + summary = paste0( + "The ", stat_name, " of all residuals is ", all$value, + if (nrow(series) > 1) paste0(", and is largest for ", worst$series, " (", worst$value, ")"), + ". ", trend + ) + ) + stats::setNames(tags, paste0("joint.residuals.", names(tags))) +} + +# Key quantities recorded in key_quantities.csv (joint.residuals.*) +joint_residuals_kqs <- function(resid, stat, scaled) { + by_index <- stat[stat$series != "All indices", , drop = FALSE] + value <- function(rows) { + if (scaled) sub("^SDNR = ", "", rows$text) else rows$value + } + kqs <- list( + start.year = min(resid$year), + end.year = max(resid$year), + n.indices = nrow(by_index), + residual.type = if (scaled) "standardized" else "log", + stat.all = stat$text[stat$series == "All indices"], + stat.by.index = paste0(by_index$series, ": ", value(by_index), collapse = "; ") + ) + stats::setNames(kqs, paste0("joint.residuals.", names(kqs))) +} diff --git a/R/plot_runs_test.R b/R/plot_runs_test.R new file mode 100644 index 00000000..8fcdf442 --- /dev/null +++ b/R/plot_runs_test.R @@ -0,0 +1,348 @@ +#' Plot residual runs tests for indices of abundance +#' +#' Plot the residuals of the fit to each index over time with three-sigma +#' control limits and a runs test for non-random residuals (Carvalho et al. +#' 2017), as in `JABBA` and `ss3diags` (`SSplotRunstest()`). +#' +#' @inheritParams plot_index +#' @param dat A tibble from [convert_output()] with labels "index_observed" +#' and "index_predicted". +#' @param fleet Character vector. Fleet(s) to plot. Default: NULL (all fleets) +#' @param ncol Number. Number of columns of panels. Default: NULL +#' +#' @returns A `ggplot2` object. +#' +#' @details Residuals are \eqn{\log(o) - \log(p)} for observed \eqn{o} and +#' predicted \eqn{p}, assuming lognormal errors. The one-sided runs test +#' checks whether residual signs are random over time (normal +#' approximation, as in `ss3diags`). The shaded band spans the three-sigma +#' limits, with sigma from the average moving range; it is green when the +#' test passes (p >= 0.05) and red when it fails. Residuals outside the +#' limits are filled red. +#' +#' No observations are dropped silently: a fleet observed in more than one +#' area, season, subseason, month, or sex becomes separate series (e.g., +#' "Survey (season 2)"), observations without a positive observed and +#' predicted value are reported in a message, and output from more than one +#' model or with values that cannot be told apart is an error. Years with no +#' observed index (BAM reports every year) are not observations. +#' +#' @references Carvalho, F., et al. 2017. Can diagnostic tests help identify +#' model misspecification in integrated stock assessments? Fisheries Research +#' 192: 28-40. +#' +#' @seealso [plot_index()], [convert_output()] +#' @export +#' +#' @examples +#' plot_runs_test(stockplotr::example_data, interactive = FALSE) +plot_runs_test <- function( + dat, + fleet = NULL, + module = NULL, + interactive = TRUE, + ncol = NULL, + make_rda = FALSE, + figures_dir = getwd() +) { + check_data(dat) + resid <- index_residuals(dat, fleet = fleet, module = module, interactive = interactive) + + ctl <- resid |> + dplyr::arrange(.data$series, .data$year) |> + dplyr::group_by(.data$series) |> + dplyr::summarise( + lcl = residual_control(.data$residual)$lcl, + ucl = residual_control(.data$residual)$ucl, + p = residual_control(.data$residual)$p, + xmin = min(.data$year) - 0.5, + xmax = max(.data$year) + 0.5, + .groups = "drop" + ) |> + dplyr::mutate( + result = factor( + dplyr::case_when( + is.na(.data$p) ~ "Too few residuals", + .data$p < 0.05 ~ "Failed (p < 0.05)", + TRUE ~ "Passed" + ), + levels = c("Passed", "Failed (p < 0.05)", "Too few residuals") + ), + panel = paste0( + wrap_label(.data$series), "\n", + dplyr::if_else( + is.na(.data$p), + "p = NA", + paste0("p = ", format.pval(round(.data$p, 3), eps = 0.001)) + ) + ) + ) + resid <- resid |> + dplyr::left_join(ctl[, c("series", "lcl", "ucl", "panel")], by = "series") |> + dplyr::mutate( + outside = !is.na(.data$lcl) & (.data$residual < .data$lcl | .data$residual > .data$ucl) + ) + + plt <- ggplot2::ggplot(resid, ggplot2::aes(x = .data$year, y = .data$residual)) + + ggplot2::geom_rect( + data = ctl, + ggplot2::aes( + xmin = .data$xmin, xmax = .data$xmax, + ymin = .data$lcl, ymax = .data$ucl, + fill = .data$result + ), + inherit.aes = FALSE, + alpha = 0.35, + na.rm = TRUE + ) + + ggplot2::geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") + + ggplot2::geom_segment(ggplot2::aes(xend = .data$year, yend = 0), color = "gray40") + + ggplot2::geom_point(fill = "white", shape = 21, size = 2, stroke = 0.8) + + ggplot2::geom_point(data = resid[resid$outside, ], color = "#B2182B", size = 2) + + ggplot2::scale_fill_manual( + values = c( + "Passed" = "#4DAC26", + "Failed (p < 0.05)" = "#D7191C", + "Too few residuals" = "gray80" + ), + drop = TRUE + ) + + ggplot2::labs(x = "Year", y = "Residual (log observed / predicted)", fill = "Runs test") + + theme_noaa() + + ggplot2::theme(legend.position = "bottom") + + ggplot2::facet_wrap(ggplot2::vars(.data$panel), ncol = ncol, scales = "free") + + if (make_rda) { + fit_export_rda( + plt, + topic_label = "runs_test", + tags = runs_test_tags(resid, ctl), + kqs = runs_test_kqs(resid, ctl), + dat = dat, + figures_dir = figures_dir + ) + } + plt +} + +#------------------------------------------------------------------------------ +# Index residuals +#------------------------------------------------------------------------------ + +# Columns that distinguish observations of a fleet in the same year +index_series_cols <- c("area", "season", "subseason", "month", "sex") + +# Observed and predicted index with log residuals, one row per series and +# year. A series is a fleet, or a fleet in one area, season, subseason, +# month, or sex when it has observations in more than one of these. +index_residuals <- function(dat, fleet = NULL, module = NULL, interactive = TRUE) { + idx <- filter_data( + dat, + label_name = "^index_(observed|predicted)$", + module = module, + geom = "point", + interactive = interactive + ) + models <- unique(stats::na.omit(as.character(idx$model))) + if (length(models) > 1) { + cli::cli_abort(c( + "{.arg dat} has output from {length(models)} models.", + "i" = "Diagnostic plots show one model at a time; supply the output of one model." + )) + } + if (!is.null(fleet)) { + idx <- idx[idx$fleet %in% fleet, , drop = FALSE] + if (nrow(idx) == 0) { + cli::cli_abort("No index found for fleet {.val {fleet}}.") + } + } + ids <- intersect(index_series_cols, names(idx)) + type <- tolower(as.character(idx$uncertainty_label)) + long <- tibble::tibble( + fleet = as.character(idx$fleet), + year = as.numeric(idx$year), + kind = sub("^index_", "", idx$label), + value = as.numeric(idx$estimate), + se = dplyr::case_when( + type %in% c("se", "sd", "log_sd", "logsd", "stddev") ~ as.numeric(idx$uncertainty), + type == "cv" ~ sqrt(log(1 + as.numeric(idx$uncertainty)^2)), + TRUE ~ NA_real_ + ) + ) + for (col in ids) long[[col]] <- as.character(idx[[col]]) + long$fleet[is.na(long$fleet)] <- "All" + # Rows that are exact copies are combined; other duplicates are an error + long <- dplyr::distinct(long) + keys <- c("fleet", "year", "kind", ids) + dups <- long |> + dplyr::group_by(dplyr::across(dplyr::all_of(keys))) |> + dplyr::filter(dplyr::n() > 1) |> + dplyr::ungroup() |> + dplyr::distinct(.data$fleet, .data$year) + if (nrow(dups) > 0) { + examples <- utils::head(paste(dups$fleet, dups$year), 3) + cli::cli_abort(c( + "{.arg dat} has more than one index value for {nrow(dups)} fleet{?s} and year{?s} that cannot be told apart (e.g., {examples}).", + "i" = "Filter {.arg dat} to one value per fleet and year, or add a column (e.g., {.field season}) that distinguishes them." + )) + } + wide <- tidyr::pivot_wider( + long, + id_cols = dplyr::all_of(c("fleet", "year", ids)), + names_from = "kind", + values_from = c("value", "se") + ) + if (!"value_predicted" %in% names(wide)) wide$value_predicted <- NA_real_ + if (!"se_observed" %in% names(wide)) wide$se_observed <- NA_real_ + out <- tibble::tibble( + fleet = wide$fleet, + year = wide$year, + observed = wide$value_observed, + predicted = wide$value_predicted, + se = wide$se_observed + ) + # Name series after the columns that vary within a fleet + extra <- rep("", nrow(out)) + for (col in ids) { + vals <- wide[[col]] + varies <- tapply(vals, out$fleet, function(v) length(unique(stats::na.omit(v))) > 1) + use <- out$fleet %in% names(varies)[varies %in% TRUE] & !is.na(vals) + extra[use] <- paste0(extra[use], ifelse(extra[use] == "", "", ", "), col, " ", vals[use]) + } + out$series <- ifelse(extra == "", out$fleet, paste0(out$fleet, " (", extra, ")")) + out <- out[!is.na(out$observed), , drop = FALSE] + out$residual <- suppressWarnings(log(out$observed) - log(out$predicted)) + out$residual[!is.finite(out$residual)] <- NA_real_ + no_resid <- is.na(out$residual) + if (any(no_resid)) { + cli::cli_alert_info( + "{sum(no_resid)} observation{?s} without a positive observed and predicted index {?has/have} no residual and {?is/are} not shown." + ) + } + out <- out[!no_resid, , drop = FALSE] + if (nrow(out) == 0) { + cli::cli_abort("No index residuals found in {.arg dat}.") + } + dplyr::arrange(out, .data$series, .data$year) +} + +# Wrap long labels for panel titles at spaces or underscores (fleet names +# such as "Survey_1_shelikof_acoustic" often have no spaces), keeping every +# character of the label +wrap_label <- function(x, width = 28) { + vapply(x, function(lab) { + if (nchar(lab) <= width) { + return(lab) + } + # Split after each space or underscore, keeping the separator + pieces <- regmatches(lab, gregexpr("[^ _]+[ _]?", lab))[[1]] + lines <- character(0) + current <- "" + for (piece in pieces) { + if (nchar(current) > 0 && nchar(current) + nchar(piece) > width) { + lines <- c(lines, trimws(current, "right")) + current <- piece + } else { + current <- paste0(current, piece) + } + } + paste(c(lines, trimws(current, "right")), collapse = "\n") + }, character(1), USE.NAMES = FALSE) +} + +#------------------------------------------------------------------------------ +# Runs test and control limits +#------------------------------------------------------------------------------ + +# Three-sigma control limits and runs test for a series of residuals +# (Carvalho et al. 2017; Winker et al. 2018). The standard deviation is +# estimated from the average moving range (Montgomery 2009), after removing +# moving ranges greater than the upper limit for moving ranges. +residual_control <- function(x) { + x <- x[!is.na(x)] + if (length(x) < 2) { + return(list(lcl = NA_real_, ucl = NA_real_, p = residual_runs_test(x))) + } + mr <- abs(diff(x)) + amr <- mean(mr) + mr <- mr[mr < 3.267 * amr] + stdev <- if (length(mr) > 0) mean(mr) / 1.128 else NA_real_ + list( + lcl = -3 * stdev, + ucl = 3 * stdev, + p = residual_runs_test(x) + ) +} + +# One-sided Wald-Wolfowitz runs test on the signs of x (in order). Returns +# the p-value for fewer runs than expected (positive autocorrelation), using +# the normal approximation. +residual_runs_test <- function(x) { + s <- sign(x[x != 0 & !is.na(x)]) + n1 <- sum(s > 0) + n2 <- sum(s < 0) + n <- n1 + n2 + if (n < 3) { + return(NA_real_) + } + # A single run: probability that all n signs are the same by chance + if (n1 == 0 || n2 == 0) { + return(2 * 0.5^n) + } + runs <- 1 + sum(s[-1] != s[-length(s)]) + mu <- 2 * n1 * n2 / n + 1 + v <- 2 * n1 * n2 * (2 * n1 * n2 - n) / (n^2 * (n - 1)) + if (v <= 0) { + return(NA_real_) + } + stats::pnorm((runs - mu) / sqrt(v)) +} + +#------------------------------------------------------------------------------ +# Caption and alternative text tags +#------------------------------------------------------------------------------ + +# `resid` has one row per residual and `ctl` one row per series +runs_test_tags <- function(resid, ctl) { + passed <- ctl$series[!is.na(ctl$p) & ctl$p >= 0.05] + failed <- ctl$series[!is.na(ctl$p) & ctl$p < 0.05] + outside <- resid[resid$outside, , drop = FALSE] + outside_text <- paste0(outside$series, " in ", outside$year) + if (length(outside_text) > 5) { + outside_text <- c(outside_text[1:5], paste(length(outside_text) - 5, "more")) + } + tags <- list( + panel.desc = fit_panel_text(ctl$series, what = "index"), + start.year = min(resid$year), + end.year = max(resid$year), + summary = paste0( + "The runs test passed for ", length(passed), " of ", nrow(ctl), " indices", + if (length(failed) > 0) paste0(" and failed for ", fit_list_text(failed)), ". ", + if (nrow(outside) == 0) { + "No residuals are outside the control limits." + } else { + paste0( + if (nrow(outside) == 1) "One residual is" else paste(nrow(outside), "residuals are"), + " outside the control limits: ", fit_list_text(outside_text), "." + ) + } + ) + ) + stats::setNames(tags, paste0("runs.test.", names(tags))) +} + +# Key quantities recorded in key_quantities.csv (runs.test.*) +runs_test_kqs <- function(resid, ctl) { + tested <- ctl[!is.na(ctl$p), , drop = FALSE] + failed <- tested$series[tested$p < 0.05] + kqs <- list( + start.year = min(resid$year), + end.year = max(resid$year), + n.indices = nrow(ctl), + n.passed = sum(tested$p >= 0.05), + failed = if (length(failed) == 0) "none" else paste(failed, collapse = "; "), + p.values = paste0(ctl$series, ": ", ifelse(is.na(ctl$p), "NA", format(round(ctl$p, 3), nsmall = 3)), collapse = "; "), + n.outside = sum(resid$outside) + ) + stats::setNames(kqs, paste0("runs.test.", names(kqs))) +} diff --git a/R/utils.R b/R/utils.R index 4688f4fe..551f8eaa 100644 --- a/R/utils.R +++ b/R/utils.R @@ -1,6 +1,6 @@ -# @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +############################################################################ # General utility functions -# @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +############################################################################ # Check end_year isn't past current year for non-projections plots # make year character if not null diff --git a/R/utils_fit_captions.R b/R/utils_fit_captions.R new file mode 100644 index 00000000..307543f5 --- /dev/null +++ b/R/utils_fit_captions.R @@ -0,0 +1,114 @@ +# Captions, alternative text, and key quantities for the fit diagnostic +# plots. Each plot has a row in captions_alt_text_template.csv whose +# placeholders start with the plot's topic label (e.g., runs.test.summary), +# and rows in key_quantity_template.csv for the quantities it records (e.g., +# runs.test.n.passed). The plot describes what it shows as tags, which +# fit_export_rda() inserts with insert_kqs(), and its results as key +# quantities, which fit_export_rda() records with export_kqs(). + +# Export key quantities, insert tags into captions_alt_text.csv, and export +# the figure as an rda +fit_export_rda <- function(plt, topic_label, tags, dat, figures_dir, kqs = list()) { + if (length(kqs) > 0) { + fit_add_kq_rows(names(kqs), dir = figures_dir) + fit_call_with_names(export_kqs, lapply(kqs, fit_tag_value), dir = figures_dir) + } + fit_add_caption_rows(topic_label, dir = figures_dir) + fit_call_with_names(insert_kqs, lapply(tags, fit_tag_value), dir = figures_dir) + create_rda( + object = plt, + topic_label = topic_label, + fig_or_table = "figure", + dat = dat, + dir = figures_dir + ) + invisible(tags) +} + +# Call export_kqs() or insert_kqs(), which take the names of the quantities +# from the names of their arguments, with one variable per value +fit_call_with_names <- function(fun, values, dir) { + env <- list2env(values, parent = environment()) + do.call(fun, c(rlang::syms(names(values)), list(dir = dir)), envir = env) +} + +# Add the template rows for key quantities to an existing key_quantities.csv +# that was made before they existed (e.g., by an older version of +# stockplotr), so that their values are recorded +fit_add_kq_rows <- function(kq_names, dir = getwd()) { + path <- fs::path(dir, "key_quantities.csv") + if (!file.exists(path)) { + return(invisible(NULL)) + } + kqs <- utils::read.csv(path) + missing <- setdiff(kq_names, kqs$key_quantity) + if (length(missing) == 0) { + return(invisible(NULL)) + } + template <- utils::read.csv( + system.file("resources", "key_quantity_template.csv", package = "stockplotr") + ) + utils::write.csv( + rbind(kqs, template[template$key_quantity %in% missing, , drop = FALSE]), + path, + row.names = FALSE + ) + invisible(NULL) +} + +# Add the template rows for a topic to an existing captions_alt_text.csv +# that was made before the topic existed (e.g., by an older version of +# stockplotr), so that its caption and alternative text are filled in +fit_add_caption_rows <- function(topic_label, dir = getwd()) { + path <- fs::path(dir, "captions_alt_text.csv") + if (!file.exists(path)) { + return(invisible(NULL)) + } + caps <- utils::read.csv(path) + if (any(caps$label == topic_label & caps$type == "figure")) { + return(invisible(NULL)) + } + template <- utils::read.csv( + system.file("resources", "captions_alt_text_template.csv", package = "stockplotr") + ) + new_rows <- template[template$label == topic_label & template$type == "figure", , drop = FALSE] + utils::write.csv(rbind(caps, new_rows), path, row.names = FALSE) + invisible(NULL) +} + +# One string per tag. Remove characters with special meaning in the +# replacement strings of stringr::str_replace_all(). +fit_tag_value <- function(x) { + x <- paste(x, collapse = " ") + gsub("[$\\\\]", "", x) +} + +#------------------------------------------------------------------------------ +# Text helpers +#------------------------------------------------------------------------------ + +# "A", "A and B", or "A, B, and C" +fit_list_text <- function(x) { + x <- unique(as.character(x)) + n <- length(x) + if (n == 0) { + return("") + } + if (n <= 2) { + return(paste(x, collapse = " and ")) + } + paste0(paste(x[-n], collapse = ", "), ", and ", x[n]) +} + +# Description of the panels in a faceted plot +fit_panel_text <- function(panels, what = "fleet") { + panels <- unique(panels) + if (length(panels) == 1) { + paste0("A single panel for ", panels, ".") + } else { + paste0( + "A set of ", length(panels), " panels, one for each ", what, + " (", fit_list_text(panels), ")." + ) + } +} diff --git a/inst/WORDLIST b/inst/WORDLIST index c725269d..b51950b8 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -3,6 +3,7 @@ BMSY Barbeaux Birthseason CPUE +Carvalho Dorn FIMS FLIMIT @@ -29,6 +30,7 @@ SPR Spawners StockSmart Subseason +al alttext asar assmt @@ -44,6 +46,7 @@ dat dataframe's diagonostic dplyr +et facetting fims flextable @@ -56,6 +59,8 @@ is’ kableExtra kobe len +loess +misspecification msy natively noaa @@ -64,11 +69,13 @@ pkgdown preformatted rceattle rda +rds rda's rdat rowid satf sd +SDNR se seealso socio diff --git a/inst/resources/captions_alt_text_template.csv b/inst/resources/captions_alt_text_template.csv index e5cfabab..e7bac6ce 100644 --- a/inst/resources/captions_alt_text_template.csv +++ b/inst/resources/captions_alt_text_template.csv @@ -48,3 +48,5 @@ catchability,table,Catchability over time for fleet catchability.fleet., sensitivity.runs,table,"Forecasted catch, biomass, spawning biomass, recruitment, fishing mortality, and stock status as projected from the base model configuration.", catch_comp,figure,Fishery age composition (caa.start.year-caa.end.year). The area of the circle is proportional to the catch. Diagonal lines indicated the top 5\\% strongest year classes.,"Bubble plot showing estimated total catch at age. The x axis shows the year, which spans from caa.start.year to caa.end.year. The y axis shows age, which spans from caa.age.min to caa.age.max. The size of the bubbles range from tot.catch.min to tot.catch.max." age_to_length_conversion,figure,"Age-to-length conversion shown as the probability distribution of length bins at each age. Density curves show the proportion of fish in each length bin at a given age, with wider curves indicating greater variability in length at age.","A set of density plots showing the age-to-length conversion for species, with one panel for each morph. The x axis shows age in years and the y axis shows length bins (in length.units). At each age, a vertical density curve shows the distribution of lengths. Mean length increases from about length.at.min.age length.units at age min.age to about length.at.max.age length.units at age max.age, after which the distributions stabilize, indicating that growth slows as fish age. " +runs_test,figure,"Residuals of the fit to indices of abundance over time, with three-sigma control limits (shaded) and a runs test for non-random residuals. The shading is green when the residuals pass the runs test (p >= 0.05) and red when they fail, and points filled in red are outside the control limits. Residuals are the log of the observed index divided by the predicted index.","runs.test.panel.desc The x axis shows year, from runs.test.start.year to runs.test.end.year, and the y axis shows the residual, with a vertical line from zero to each residual. runs.test.summary" +joint_residuals,figure,"joint.residuals.type of the fits to joint.residuals.index.list joint.residuals.layout.","The x axis shows year, from joint.residuals.start.year to joint.residuals.end.year, and the y axis shows the residual, with a vertical line from zero to each residual, colored by index. joint.residuals.summary" diff --git a/inst/resources/fims_var_names.csv b/inst/resources/fims_var_names.csv index 3d295826..7b3c2cd1 100644 --- a/inst/resources/fims_var_names.csv +++ b/inst/resources/fims_var_names.csv @@ -1,46 +1,59 @@ -module_name,label,alt_label -Fleet,age_to_length_conversion, -Fleet,agecomp_expected,age_composition_expected -Fleet,agecomp_proportion,age_composition_proportion -Fleet,index_expected, -Fleet,index_naa,index -Fleet,index_nal,index -Fleet,index_numbers, -Fleet,index_waa,index_weight -Fleet,index_weight, -Fleet,landings_expected, -Fleet,landings_naa,landings -Fleet,landings_nal,landings -Fleet,landings_numbers,landings_numbers -Fleet,landings_waa,landings_weight -Fleet,landings_weight,landings_weight -Fleet,lengthcomp_expected,length_composition_expected -Fleet,lengthcomp_proportion,length_composition_proportion -Fleet,log_Fmort,log_fishing_mortality -Fleet,log_q,log_catchability -Population,Biomass,biomass -Population,NAA,abundance -Population,Recruitment,recruitment -Population,SSB,spawning_biomass -Population,log_M,log_natural_mortality -Population,log_init_naa,log_abundance_init -maturity,inflection_point, -maturity,slope, -recruitment,log_devs,log_recruitment_deviations -recruitment,log_rzero,log_recruitment_unfished -recruitment,logit_steep, -selectivity,inflection_point, -selectivity,slope, -,AgeCompositionExpected,age_composition_expected -,AgeCompositionProportion,age_composition_proportion -,FMort,fishing_mortality -,IndexExpected,index_expected -,IndexNumberAtAge,index_numbers -,IndexNumberAtLength,index_numbers -,LandingsExpected,landings_expected -,LandingsNumberAtAge,landings_numbers -,LandingsNumberAtLength,landings_numbers -,LengthCompositionExpected,length_composition_expected -,LengthCompositionProportion,length_composition_proportion -,LogRecDev,log_recruitment_deviations -,Q,catchability +module_name,label,alt_label,unit +Fleet,age_to_length_conversion,, +Fleet,agecomp_expected,age_composition_expected, +Fleet,agecomp_proportion,age_composition_proportion, +Fleet,index_expected,, +Fleet,index_naa,index,numbers +Fleet,index_nal,index,numbers +Fleet,index_numbers,,numbers +Fleet,index_waa,index_weight,mt +Fleet,index_weight,,mt +Fleet,landings_expected,, +Fleet,landings_naa,landings,numbers +Fleet,landings_nal,landings,numbers +Fleet,landings_numbers,landings_numbers,numbers +Fleet,landings_waa,landings_weight,mt +Fleet,landings_weight,landings_weight,mt +Fleet,lengthcomp_expected,length_composition_expected, +Fleet,lengthcomp_proportion,length_composition_proportion, +Fleet,log_Fmort,log_fishing_mortality, +Fleet,log_q,log_catchability, +Population,Biomass,biomass,mt +Population,NAA,abundance,numbers +Population,Recruitment,recruitment,numbers +Population,SSB,spawning_biomass,mt +Population,log_M,log_natural_mortality, +Population,log_init_naa,log_abundance_init, +maturity,inflection_point,, +maturity,slope,, +recruitment,log_devs,log_recruitment_deviations, +recruitment,log_rzero,log_recruitment_unfished, +recruitment,logit_steep,, +selectivity,inflection_point,, +selectivity,slope,, +,AgeCompositionExpected,age_composition_expected, +,AgeCompositionProportion,age_composition_proportion, +,FMort,fishing_mortality, +,IndexExpected,index_expected, +,IndexNumberAtAge,index_numbers, +,IndexNumberAtLength,index_numbers, +,LandingsExpected,landings_expected, +,LandingsNumberAtAge,landings_numbers, +,LandingsNumberAtLength,landings_numbers, +,LengthCompositionExpected,length_composition_expected, +,LengthCompositionProportion,length_composition_proportion, +,LogRecDev,log_recruitment_deviations, +,Q,catchability, +Population,expected_recruitment,recruitment,numbers +Population,numbers_at_age,abundance,numbers +Population,unfished_spawning_biomass,spawning_biomass_unfished,mt +Population,unfished_biomass,biomass_unfished,mt +Population,mortality_M,natural_mortality, +Recruitment,log_devs,log_recruitment_deviations, +Recruitment,log_rzero,log_recruitment_unfished, +Recruitment,log_sd,log_recruitment_sd, +Recruitment,logit_steep,logit_steepness, +Population,biomass,,mt +Population,spawning_biomass,,mt +Fleet,landings_observed,,mt +Fleet,landings_predicted,,mt diff --git a/inst/resources/key_quantity_template.csv b/inst/resources/key_quantity_template.csv index 04fb14ae..d2ecfc46 100644 --- a/inst/resources/key_quantity_template.csv +++ b/inst/resources/key_quantity_template.csv @@ -53,6 +53,12 @@ fecundity.units,,units for fecundity (number of eggs),,,, fork.length.max,,maximum fork length,,,, fork.length.min,,minimum fork length,,,, fork.length.units,,units for fork length,,,, +joint.residuals.end.year,,final year of the index residuals in the joint residuals plot,,,, +joint.residuals.n.indices,,number of indices in the joint residuals plot,,,, +joint.residuals.residual.type,,type of residuals in the joint residuals plot (log or standardized),,,, +joint.residuals.start.year,,start year of the index residuals in the joint residuals plot,,,, +joint.residuals.stat.all,,fit statistic of all index residuals (RMSE of log residuals or SDNR of standardized residuals with its critical value),,,, +joint.residuals.stat.by.index,,fit statistic of the residuals of each index,,,, kobe.end.year,,final year displayed on the Kobe plot,,,, landings.end.year,,final year of historical landings data,,,, landings.max,,maximum historical landings,,,, @@ -112,6 +118,13 @@ rel.recruitment.max,,maximum relative recruitment (R/R0),,,, rel.recruitment.min,,minimum relative recruitment (R/R0),,,, rel.ssb.max,,maximum relative spawning biomass,,,, rel.ssb.min,,minimum relative spawning biomass,,,, +runs.test.end.year,,final year of the index residuals in the runs test plot,,,, +runs.test.failed,,indices whose residuals fail the runs test (p < 0.05),,,, +runs.test.n.indices,,number of indices in the runs test plot,,,, +runs.test.n.outside,,number of index residuals outside the three-sigma control limits,,,, +runs.test.n.passed,,number of indices whose residuals pass the runs test (p >= 0.05),,,, +runs.test.p.values,,p-value of the runs test for each index,,,, +runs.test.start.year,,start year of the index residuals in the runs test plot,,,, selectivity.end.year,,final year for selectivity estimates,,,, selectivity.start.year,,start year for selectivity estimates,,,, spr.start.year,,start year of SPR timeseries,,,, diff --git a/man/conout_bam.Rd b/man/conout_bam.Rd index 7b3ee972..8db1ca02 100644 --- a/man/conout_bam.Rd +++ b/man/conout_bam.Rd @@ -8,7 +8,8 @@ conout_bam(file, out_new, fleet_names = NULL) } \arguments{ \item{file}{Path. Assessment model output file path. Currently accepts SS3 -\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object.} +\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output +(the tibble returned by \code{FIMS::get_estimates()}, a \code{FIMSFit} object, or an .rds file with either).} \item{fleet_names}{Character vector. Names of fleets in the assessment model as shortened in the output file or names of fleets to replace the native fleet diff --git a/man/conout_rceattle.Rd b/man/conout_rceattle.Rd index 85612b51..e053b2a8 100644 --- a/man/conout_rceattle.Rd +++ b/man/conout_rceattle.Rd @@ -8,7 +8,8 @@ conout_rceattle(file, out_new, fleet_names = NULL) } \arguments{ \item{file}{Path. Assessment model output file path. Currently accepts SS3 -\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object.} +\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output +(the tibble returned by \code{FIMS::get_estimates()}, a \code{FIMSFit} object, or an .rds file with either).} \item{fleet_names}{Character vector. Names of fleets in the assessment model as shortened in the output file or names of fleets to replace the native fleet diff --git a/man/conout_ss3.Rd b/man/conout_ss3.Rd index fe5de1e8..303b7fbf 100644 --- a/man/conout_ss3.Rd +++ b/man/conout_ss3.Rd @@ -8,7 +8,8 @@ conout_ss3(file, out_new, fleet_names = NULL) } \arguments{ \item{file}{Path. Assessment model output file path. Currently accepts SS3 -\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object.} +\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output +(the tibble returned by \code{FIMS::get_estimates()}, a \code{FIMSFit} object, or an .rds file with either).} \item{fleet_names}{Character vector. Names of fleets in the assessment model as shortened in the output file or names of fleets to replace the native fleet diff --git a/man/convert_output.Rd b/man/convert_output.Rd index 782b50e7..b409306e 100644 --- a/man/convert_output.Rd +++ b/man/convert_output.Rd @@ -4,11 +4,18 @@ \alias{convert_output} \title{Convert Output} \usage{ -convert_output(file, model = NULL, fleet_names = NULL, save_dir = NULL) +convert_output( + file, + model = NULL, + fleet_names = NULL, + save_dir = NULL, + start_year = NULL +) } \arguments{ \item{file}{Path. Assessment model output file path. Currently accepts SS3 -\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object.} +\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output +(the tibble returned by \code{FIMS::get_estimates()}, a \code{FIMSFit} object, or an .rds file with either).} \item{model}{String. Assessment model used in evaluation @@ -28,6 +35,12 @@ saved file name is not indicated, the file will be saved as "std_output.rda" in the file path provided. Default: NULL} + +\item{start_year}{Number. First year of the model, for FIMS output only. +FIMS reports years as indices (1 = first model year), which are converted +to years when \code{start_year} is given. + +Default: NULL (years are left as indices)} } \value{ A reformatted and standardized tibble of assessment model results @@ -41,6 +54,14 @@ Format stock assessment output file to a standardized tibble. The resulting object is a long-format, transformed and machine readable version of a model output file. Converted data frame is always returned. It will also be saved if save_dir is not NULL. + +For FIMS, each data point (index, landings, and age and length +compositions) is split into an observed row (e.g., "index_observed") and a +predicted row (e.g., "index_predicted"), as for the other models. Observed +rows hold the observation error (a log-scale standard error for lognormal +data) in \code{uncertainty} and the negative log-likelihood in \code{likelihood}. +Other estimates are relabeled with \code{fims_var_names.csv} (e.g., +"expected_recruitment" becomes "recruitment"). } \examples{ \dontrun{ diff --git a/man/plot_joint_residuals.Rd b/man/plot_joint_residuals.Rd new file mode 100644 index 00000000..4dd03d44 --- /dev/null +++ b/man/plot_joint_residuals.Rd @@ -0,0 +1,120 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_joint_residuals.R +\name{plot_joint_residuals} +\alias{plot_joint_residuals} +\title{Plot joint residuals of the fits to indices of abundance} +\usage{ +plot_joint_residuals( + dat, + fleet = NULL, + module = NULL, + scaled = FALSE, + overlay = TRUE, + show_stat = TRUE, + ncol = NULL, + interactive = TRUE, + make_rda = FALSE, + figures_dir = getwd() +) +} +\arguments{ +\item{dat}{A tibble from \code{\link[=convert_output]{convert_output()}} with labels "index_observed" +and "index_predicted".} + +\item{fleet}{Character vector. Fleet(s) to plot. Default: NULL (all fleets)} + +\item{module}{Character vector. (Optional) Module name found in \code{dat$module_name}. +If selecting >1 module, place them in a vector like c("module1", "module2"). + +Default: NULL + +If the interactive and >1 module_name is found, user will select the +module_name in the console. @seealso \code{\link[=filter_data]{filter_data()}}} + +\item{scaled}{Logical. TRUE/FALSE; plot standardized residuals (log +residuals divided by their standard error) instead of log residuals. + +Default: FALSE} + +\item{overlay}{Logical. TRUE/FALSE; plot all indices on one panel, with a +boxplot of the residuals in each year. FALSE plots each index on its own +panel. + +Default: TRUE} + +\item{show_stat}{Logical. TRUE/FALSE; show the fit statistic in the top +right corner of each panel: the root mean square error (RMSE) of log +residuals, or the standard deviation of normalized residuals (SDNR) of +scaled residuals. When \code{overlay = TRUE}, the statistic of each index is +also shown in the legend. + +Default: TRUE} + +\item{ncol}{Number. Number of columns of panels when \code{overlay = FALSE}. + +Default: NULL} + +\item{interactive}{Logical. TRUE/FALSE; indicate whether the environment is interactive. + +Default: \code{FALSE}} + +\item{make_rda}{Logical. TRUE/FALSE; indicate whether to save the object and +make an automated caption and alternative text in the form of an \code{rda} object. If TRUE, +the rda will be exported to the folder indicated in the argument "figures_dir". + +Default: \code{FALSE}.} + +\item{figures_dir}{Path. Path to the "figures" folder. + +Default: \code{getwd()} + +The folder is created within the path if it does not exist.} +} +\value{ +A \code{ggplot2} object. +} +\description{ +Plot the residuals of the fits to all indices over time with a loess +smoother and a fit statistic, as in \code{JABBA} (\code{jbplot_residuals()} and +\code{jbplot_stdresiduals()}) and \code{ss3diags} (\code{SSplotJABBAres()}). A trend in +the smoother shows a pattern in the residuals; a large statistic shows a +poor fit. +} +\details{ +Log residuals are \eqn{r = \log(o) - \log(p)} for observed +\eqn{o} and predicted \eqn{p}, as in \code{\link[=plot_runs_test]{plot_runs_test()}}; scaled +residuals are \eqn{z = r / s} for the log-scale standard error \eqn{s} of +the observation. The RMSE is \eqn{100 \sqrt{\overline{r^2}}} percent. The +SDNR is \eqn{\sqrt{\sum z^2 / (n - 1)}}, shown with its critical value +\eqn{\sqrt{\chi^2_{0.95, n - 1} / (n - 1)}}; an SDNR above the critical +value suggests the standard errors are too small (Francis 2011). + +The smoother is a loess regression (span 0.75) of the mean residual of +each year on year, for all indices together or for each index on its own +panel. It is not drawn for fewer than five years of residuals, or when +the years are too sparse for a stable fit; a message names the indices +without one. + +Observations are treated as in \code{\link[=plot_runs_test]{plot_runs_test()}}: a fleet observed in +more than one area, season, subseason, month, or sex becomes separate +series, and observations without a positive observed and predicted value +are reported in a message. Scaled residuals also need a standard error; +observations without one are reported in a message. +} +\examples{ +plot_joint_residuals(stockplotr::example_data, interactive = FALSE) +plot_joint_residuals( + stockplotr::example_data, + scaled = TRUE, + overlay = FALSE, + interactive = FALSE +) +} +\references{ +Francis, R.I.C.C. 2011. Data weighting in statistical +fisheries stock assessment models. Canadian Journal of Fisheries and +Aquatic Sciences 68: 1124-1138. +} +\seealso{ +\code{\link[=plot_runs_test]{plot_runs_test()}}, \code{\link[=plot_index]{plot_index()}}, \code{\link[=convert_output]{convert_output()}} +} diff --git a/man/plot_runs_test.Rd b/man/plot_runs_test.Rd new file mode 100644 index 00000000..a3687cac --- /dev/null +++ b/man/plot_runs_test.Rd @@ -0,0 +1,83 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_runs_test.R +\name{plot_runs_test} +\alias{plot_runs_test} +\title{Plot residual runs tests for indices of abundance} +\usage{ +plot_runs_test( + dat, + fleet = NULL, + module = NULL, + interactive = TRUE, + ncol = NULL, + make_rda = FALSE, + figures_dir = getwd() +) +} +\arguments{ +\item{dat}{A tibble from \code{\link[=convert_output]{convert_output()}} with labels "index_observed" +and "index_predicted".} + +\item{fleet}{Character vector. Fleet(s) to plot. Default: NULL (all fleets)} + +\item{module}{Character vector. (Optional) Module name found in \code{dat$module_name}. +If selecting >1 module, place them in a vector like c("module1", "module2"). + +Default: NULL + +If the interactive and >1 module_name is found, user will select the +module_name in the console. @seealso \code{\link[=filter_data]{filter_data()}}} + +\item{interactive}{Logical. TRUE/FALSE; indicate whether the environment is interactive. + +Default: \code{FALSE}} + +\item{ncol}{Number. Number of columns of panels. Default: NULL} + +\item{make_rda}{Logical. TRUE/FALSE; indicate whether to save the object and +make an automated caption and alternative text in the form of an \code{rda} object. If TRUE, +the rda will be exported to the folder indicated in the argument "figures_dir". + +Default: \code{FALSE}.} + +\item{figures_dir}{Path. Path to the "figures" folder. + +Default: \code{getwd()} + +The folder is created within the path if it does not exist.} +} +\value{ +A \code{ggplot2} object. +} +\description{ +Plot the residuals of the fit to each index over time with three-sigma +control limits and a runs test for non-random residuals (Carvalho et al. +2017), as in \code{JABBA} and \code{ss3diags} (\code{SSplotRunstest()}). +} +\details{ +Residuals are \eqn{\log(o) - \log(p)} for observed \eqn{o} and +predicted \eqn{p}, assuming lognormal errors. The one-sided runs test +checks whether residual signs are random over time (normal +approximation, as in \code{ss3diags}). The shaded band spans the three-sigma +limits, with sigma from the average moving range; it is green when the +test passes (p >= 0.05) and red when it fails. Residuals outside the +limits are filled red. + +No observations are dropped silently: a fleet observed in more than one +area, season, subseason, month, or sex becomes separate series (e.g., +"Survey (season 2)"), observations without a positive observed and +predicted value are reported in a message, and output from more than one +model or with values that cannot be told apart is an error. Years with no +observed index (BAM reports every year) are not observations. +} +\examples{ +plot_runs_test(stockplotr::example_data, interactive = FALSE) +} +\references{ +Carvalho, F., et al. 2017. Can diagnostic tests help identify +model misspecification in integrated stock assessments? Fisheries Research +192: 28-40. +} +\seealso{ +\code{\link[=plot_index]{plot_index()}}, \code{\link[=convert_output]{convert_output()}} +} diff --git a/tests/testthat/fixtures/fims_output/README.md b/tests/testthat/fixtures/fims_output/README.md new file mode 100644 index 00000000..a16f69d9 --- /dev/null +++ b/tests/testthat/fixtures/fims_output/README.md @@ -0,0 +1,7 @@ +# FIMS fixture + +`fims_estimates.rds` is a subset of `FIMS::estimates_big_optimized` (the output of +`FIMS::get_estimates()` for the `data_big` example model shipped with FIMS, +NOAA-FIMS/FIMS commit 506b290). Rows kept: every data row (non-missing +`observed` or `likelihood`), all fixed and random effects, and the spawning +biomass, recruitment, biomass, index, and catch derived quantities. diff --git a/tests/testthat/fixtures/fims_output/fims_estimates.rds b/tests/testthat/fixtures/fims_output/fims_estimates.rds new file mode 100644 index 00000000..e2f03011 Binary files /dev/null and b/tests/testthat/fixtures/fims_output/fims_estimates.rds differ diff --git a/tests/testthat/helper-fit_captions.R b/tests/testthat/helper-fit_captions.R new file mode 100644 index 00000000..8d6d8306 --- /dev/null +++ b/tests/testthat/helper-fit_captions.R @@ -0,0 +1,38 @@ +# Evaluate plotting code with `make_rda = TRUE` in a temporary directory and +# return the saved rda (a list with the figure, caption, and alt text), with +# the key quantities recorded in key_quantities.csv as `kqs` (a named +# vector of the filled-in values). +# `captions` is an optional captions_alt_text.csv to start from, and `dir` the +# figures_dir the plot was saved to (default: the temporary directory). +fit_rda <- function(code, topic_label, captions = NULL, dir = NULL) { + wd <- tempfile("fit_rda") + dir.create(wd) + old <- setwd(wd) + on.exit( + { + setwd(old) + unlink(wd, recursive = TRUE) + }, + add = TRUE + ) + if (!is.null(captions)) { + utils::write.csv(captions, "captions_alt_text.csv", row.names = FALSE) + } + suppressMessages(invisible(utils::capture.output(force(code)))) + rda <- NULL + out <- if (is.null(dir)) wd else dir + load(file.path(out, "figures", paste0(topic_label, "_figure.rda"))) + kqs <- utils::read.csv(file.path(out, "key_quantities.csv")) + kqs <- kqs[!is.na(kqs$value) & kqs$value != "", , drop = FALSE] + rda$kqs <- stats::setNames(as.character(kqs$value), kqs$key_quantity) + rda +} + +# Expect that every placeholder starting with `prefix` (e.g., "runs.test.") +# was replaced in the caption and alt text +expect_tags_filled <- function(rda, prefix) { + expect_s3_class(rda$figure, "ggplot") + expect_no_match(c(rda$caption, rda$alt_text), prefix, fixed = TRUE) + expect_gt(nchar(rda$caption), 50) + expect_gt(nchar(rda$alt_text), 100) +} diff --git a/tests/testthat/test-convert_output_fims.R b/tests/testthat/test-convert_output_fims.R new file mode 100644 index 00000000..1b25ecc7 --- /dev/null +++ b/tests/testthat/test-convert_output_fims.R @@ -0,0 +1,143 @@ +# FIMS::get_estimates() output for the FIMS data_big example +# (see fixtures/fims_output/README.md) +fims_estimates <- function() { + readRDS(test_path("fixtures", "fims_output", "fims_estimates.rds")) +} +convert_fims <- function(dat = fims_estimates(), ...) { + suppressMessages(convert_output(dat, model = "fims", ...)) +} + +test_that("FIMS data points are split into observed and predicted rows", { + est <- fims_estimates() + out <- convert_fims(est, start_year = 2001) + counts <- table(out$label[grepl("_(observed|predicted)$", out$label)]) + expect_equal( + as.vector(counts[c("index_observed", "index_predicted", "landings_observed", "landings_predicted")]), + c(30, 30, 30, 30) + ) + expect_equal( + as.vector(counts[c("age_composition_observed", "length_composition_observed")]), + c(sum(est$label == "agecomp_expected" & !is.na(est$observed)), sum(est$label == "lengthcomp_expected" & !is.na(est$observed))) + ) + + rows <- est[est$label == "index_expected" & !is.na(est$observed), ] + obs <- out[out$label == "index_observed", ] + pred <- out[out$label == "index_predicted", ] + expect_equal(obs$estimate, rows$observed) + expect_equal(pred$estimate, rows$estimated) + # Observation error and negative log-likelihood of each observation + # Lognormal observation error is reported as a CV + expect_equal(obs$uncertainty, sqrt(exp(exp(rows$log_sd)^2) - 1)) + expect_equal(unique(obs$uncertainty_label), "cv") + expect_equal(obs$likelihood, -rows$likelihood) + expect_equal(unique(obs$fleet), "survey1") + expect_equal(obs$year, rows$year_i + 2000) + + # Age compositions keep their bins and counts + comp <- out[out$label == "age_composition_observed", ] + expect_true(all(!is.na(comp$age))) + expect_equal(sum(comp$estimate), sum(est$observed[est$label == "agecomp_expected"], na.rm = TRUE)) + expect_true(all(!is.na(out$length_bins[out$label == "length_composition_observed"]))) +}) + +test_that("FIMS expected values for fleets without data keep their FIMS label", { + out <- convert_fims() + # fleet1 is the fishery, which has no index + expect_false(any(out$label == "index_predicted" & out$fleet == "fleet1")) + expect_true(any(out$label == "index_expected" & out$fleet == "fleet1")) +}) + +test_that("FIMS observations coded as missing are not observations", { + est <- fims_estimates() + idx <- which(est$label == "index_expected" & !is.na(est$observed)) + est$observed[idx[1:2]] <- -999 + out <- convert_fims(est) + expect_equal(sum(out$label == "index_observed"), length(idx) - 2) + expect_equal(sum(out$label == "index_predicted"), length(idx)) +}) + +test_that("FIMS lognormal intervals are on the original scale", { + est <- fims_estimates() + est <- est[est$label == "landings_expected" | est$label == "catch_expected", ] + row <- which(!is.na(est$observed))[1] + est$observed[row] <- 100 + est$log_sd[row] <- log(0.2) + out <- convert_fims(est) + obs <- out[out$label == "landings_observed", ][1, ] + expect_equal(obs$uncertainty_label, "cv") + # The CV gives back the lognormal sigma, and so the lognormal interval + sigma <- sqrt(log(1 + obs$uncertainty^2)) + expect_equal(sigma, 0.2) + expect_equal(100 * exp(c(-1.96, 1.96) * sigma), c(67.57, 147.99), tolerance = 1e-4) +}) + +test_that("FIMS estimates that could not be calculated are NA", { + est <- fims_estimates() + ssb <- which(est$label == "spawning_biomass") + est$estimated[ssb[1]] <- -999 + est$uncertainty[ssb[1]] <- -999 + out <- convert_fims(est) + out_ssb <- out[out$label == "spawning_biomass", ] + # The row is kept, with a missing estimate + expect_equal(nrow(out_ssb), length(ssb)) + expect_true(is.na(out_ssb$estimate[1]) && is.na(out_ssb$uncertainty[1])) + expect_false(any(out$estimate %in% -999)) +}) + +test_that("FIMS estimates are relabeled and get years and eras", { + est <- fims_estimates() + out <- convert_fims(est, start_year = 2001) + ssb <- out[out$label == "spawning_biomass", ] + expect_equal(ssb$estimate, est$estimated[est$label == "spawning_biomass"]) + expect_equal(ssb$uncertainty, est$uncertainty[est$label == "spawning_biomass"]) + expect_equal(range(ssb$year), c(2001, 2031)) + # Spawning biomass at the start of the year after the model ends + expect_equal(ssb$era[ssb$year == 2031], "fore") + expect_true(all(ssb$era[ssb$year <= 2030] == "time")) + + expect_equal(sum(out$label == "recruitment"), sum(est$label == "expected_recruitment")) + expect_equal(sum(out$label == "log_recruitment_deviations"), sum(est$label == "log_devs")) + expect_false(any(c("expected_recruitment", "log_devs") %in% out$label)) + + # Without start_year, years are the FIMS year indices + expect_message(convert_output(est, model = "fims"), "start_year") + expect_equal(range(convert_fims(est)$year[convert_fims(est)$label == "spawning_biomass"]), c(1, 31)) +}) + +test_that("FIMS output is read from an rds file", { + path <- tempfile(fileext = ".rds") + saveRDS(fims_estimates(), path) + on.exit(unlink(path), add = TRUE) + expect_equal(nrow(convert_fims(path)), nrow(convert_fims())) +}) + +test_that("a FIMSFit object needs FIMS", { + skip_if(requireNamespace("FIMS", quietly = TRUE)) + methods::setClass("FIMSFit", representation(x = "numeric"), where = environment()) + expect_error( + suppressMessages(convert_output(methods::new("FIMSFit", x = 1))), + "FIMS" + ) +}) + +test_that("plots work with converted FIMS output", { + out <- convert_fims(start_year = 2001) + p <- suppressMessages(plot_index(out, interactive = FALSE)) + b <- ggplot2::ggplot_build(p) + expect_equal(nrow(b$data[[1]]), 30) + expect_equal(range(b$data[[1]]$x), c(2001, 2030)) + for (f in list(plot_landings, plot_spawning_biomass, plot_biomass, plot_recruitment, plot_recruitment_deviations)) { + p <- suppressMessages(suppressWarnings(f(out, interactive = FALSE))) + expect_s3_class(p, "ggplot") + expect_gt(sum(sapply(ggplot2::ggplot_build(p)$data, nrow)), 0) + } +}) + +test_that("plot_index plots only observed and predicted indices", { + # pollock (Rceattle) has other labels that contain "index", such as + # index_log_q, and ss3_simple has a single index fleet + for (dat in list(stockplotr::pollock, stockplotr::ss3_simple, stockplotr::example_data)) { + p <- suppressMessages(plot_index(dat, interactive = FALSE)) + expect_no_error(ggplot2::ggplot_build(p)) + } +}) diff --git a/tests/testthat/test-plot_joint_residuals.R b/tests/testthat/test-plot_joint_residuals.R new file mode 100644 index 00000000..64e1a1fb --- /dev/null +++ b/tests/testthat/test-plot_joint_residuals.R @@ -0,0 +1,95 @@ +test_that("plot_joint_residuals works with each example model output and option", { + est <- readRDS(test_path("fixtures", "fims_output", "fims_estimates.rds")) + datasets <- list( + example_data = stockplotr::example_data, + pollock = stockplotr::pollock, + ss3_simple = stockplotr::ss3_simple, + vermilion_snapper = stockplotr::vermilion_snapper, + fims = suppressMessages(convert_output(est, model = "fims", start_year = 2001)) + ) + for (name in names(datasets)) { + # ss3_simple also has catch labeled as an index in the "catch" module + module <- if (name == "ss3_simple") "cpue" else NULL + for (scaled in c(FALSE, TRUE)) { + p <- suppressMessages(plot_joint_residuals( + datasets[[name]], + module = module, scaled = scaled, overlay = !scaled, interactive = FALSE + )) + expect_no_error(ggplot2::ggplot_build(p)) + } + } +}) + +test_that("residual_stat gives the RMSE and the SDNR with its critical value", { + resid <- data.frame(series = rep(c("A", "B"), c(3, 1)), residual = c(0.1, -0.2, 0.3, 0.4)) + rmse <- residual_stat(resid) + expect_equal(rmse$series, c("A", "B", "All indices")) + expect_equal(rmse$stat, 100 * sqrt(c(mean(c(0.01, 0.04, 0.09)), 0.16, mean(resid$residual^2)))) + expect_equal(rmse$text[2], "RMSE = 40.0%") + sdnr <- residual_stat(resid, scaled = TRUE) + expect_equal(sdnr$stat[1], sqrt(0.14 / 2)) + expect_equal(sdnr$critical[1], sqrt(stats::qchisq(0.95, df = 2) / 2)) + # One residual has no SDNR + expect_equal(sdnr$text[2], "SDNR = NA") +}) + +test_that("plot_joint_residuals overlays indices with the statistic, boxplots, and smoother", { + resid <- index_residuals(stockplotr::example_data, interactive = FALSE) + stat <- residual_stat(resid) + p <- plot_joint_residuals(stockplotr::example_data, interactive = FALSE) + layer <- function(p, geom) p$layers[[which(sapply(p$layers, function(l) inherits(l$geom, geom)))]]$data + expect_setequal(levels(p$data$index), paste0(stat$series[1:2], " (", stat$value[1:2], ")")) + expect_equal(layer(p, "GeomLabel")$text, stat$text[3]) + # Boxplots only for years with more than one residual + expect_true(all(table(resid$year)[as.character(layer(p, "GeomBoxplot")$year)] > 1)) + smooth <- layer(p, "GeomLine") + expect_equal(smooth$mean, as.numeric(tapply(resid$residual, resid$year, mean)[as.character(smooth$year)])) + + p <- plot_joint_residuals(stockplotr::example_data, show_stat = FALSE, interactive = FALSE) + expect_false(any(sapply(p$layers, function(l) inherits(l$geom, "GeomLabel")))) + expect_setequal(levels(p$data$index), stat$series[1:2]) + + short <- stockplotr::example_data[stockplotr::example_data$year %in% 1995:1998, ] + expect_message(p <- plot_joint_residuals(short, interactive = FALSE), "Too few years") + expect_false(any(sapply(p$layers, function(l) inherits(l$geom, "GeomLine")))) +}) + +test_that("plot_joint_residuals scales residuals and plots indices on panels", { + resid <- index_residuals(stockplotr::example_data, interactive = FALSE) + resid <- resid[order(resid$series, resid$year), ] + p <- plot_joint_residuals(stockplotr::example_data, scaled = TRUE, overlay = FALSE, interactive = FALSE) + plotted <- p$data[order(p$data$series, p$data$year), ] + expect_equal(plotted$residual, resid$residual / resid$se) + expect_s3_class(p$facet, "FacetWrap") + expect_false(any(sapply(p$layers, function(l) inherits(l$geom, "GeomBoxplot")))) + # One statistic per panel + labels <- p$layers[[which(sapply(p$layers, function(l) inherits(l$geom, "GeomLabel")))]]$data + expect_equal(labels$text, residual_stat(plotted, scaled = TRUE)$text[1:2]) + + # Observations without a standard error are reported, not dropped silently + dat <- stockplotr::example_data + dat$uncertainty[which(dat$label == "index_observed")[1]] <- NA + expect_message( + p <- plot_joint_residuals(dat, scaled = TRUE, interactive = FALSE), + "1 residual without a standard error" + ) + expect_equal(nrow(p$data), nrow(resid) - 1) +}) + +test_that("plot_joint_residuals saves a caption and alt text", { + rda <- fit_rda( + plot_joint_residuals(stockplotr::example_data, scaled = TRUE, interactive = FALSE, make_rda = TRUE), + "joint_residuals" + ) + expect_tags_filled(rda, "joint.residuals.") + expect_match(rda$caption, "Standardized residuals of the fits to Triennial", fixed = TRUE) + expect_match(rda$alt_text, "(SDNR) of all residuals is", fixed = TRUE) + # The results are recorded in key_quantities.csv + kqs <- rda$kqs[grepl("^joint[.]residuals[.]", names(rda$kqs))] + stat <- residual_stat(rda$figure$data, scaled = TRUE) + expect_equal(kqs[["joint.residuals.n.indices"]], "2") + expect_equal(kqs[["joint.residuals.residual.type"]], "standardized") + expect_equal(kqs[["joint.residuals.stat.all"]], stat$text[3]) + expect_match(kqs[["joint.residuals.stat.by.index"]], "^[^;]+: [0-9.]+ \\(critical [0-9.]+\\); [^;]+: ") + expect_length(kqs, 6) +}) diff --git a/tests/testthat/test-plot_runs_test.R b/tests/testthat/test-plot_runs_test.R new file mode 100644 index 00000000..47e52745 --- /dev/null +++ b/tests/testthat/test-plot_runs_test.R @@ -0,0 +1,127 @@ +fims_converted <- function() { + est <- readRDS(test_path("fixtures", "fims_output", "fims_estimates.rds")) + suppressMessages(convert_output(est, model = "fims", start_year = 2001)) +} + +test_that("residual_control estimates limits and the runs test", { + ctl <- residual_control(c(1, -1, 1, -1, 1, -1)) + # All moving ranges are 2, so sigma = 2 / 1.128 + expect_equal(ctl$ucl, 3 * 2 / 1.128) + expect_equal(ctl$lcl, -ctl$ucl) + # Alternating signs have the most runs; long runs of one sign the fewest + expect_gt(ctl$p, 0.9) + expect_lt(residual_control(c(-3:-1, 1:3))$p, 0.05) + expect_true(is.na(residual_control(1)$lcl)) + expect_equal(residual_runs_test(rep(1, 5)), 2 * 0.5^5) + # Normal approximation: 6 positive, 4 negative, 5 runs + x <- c(1, 1, -1, 1, -1, -1, -1, 1, 1, 1) + mu <- 2 * 6 * 4 / 10 + 1 + v <- 2 * 6 * 4 * (2 * 6 * 4 - 10) / (10^2 * 9) + expect_equal(residual_runs_test(x), stats::pnorm((5 - mu) / sqrt(v))) +}) + +test_that("index residuals match each example model output", { + datasets <- list( + example_data = stockplotr::example_data, + pollock = stockplotr::pollock, + ss3_simple = stockplotr::ss3_simple, + vermilion_snapper = stockplotr::vermilion_snapper, + fims = fims_converted() + ) + for (name in names(datasets)) { + dat <- datasets[[name]] + # ss3_simple also has catch labeled as an index in the "catch" module + module <- if (name == "ss3_simple") "cpue" else NULL + res <- suppressMessages(index_residuals(dat, module = module, interactive = FALSE)) + rows <- dat[dat$label %in% c("index_observed", "index_predicted"), ] + if (!is.null(module)) rows <- rows[rows$module_name == module, ] + obs <- rows[rows$label == "index_observed" & !is.na(rows$estimate), ] + pred <- rows[rows$label == "index_predicted", ] + obs <- obs[order(obs$fleet, obs$year), ] + pred <- pred[match(paste(obs$fleet, obs$year), paste(pred$fleet, pred$year)), ] + expect_equal(nrow(res), nrow(obs), label = name) + expect_equal(res$residual, log(obs$estimate) - log(pred$estimate), label = name) + # Log-scale standard errors, with CVs (BAM and FIMS) converted + se <- ifelse(obs$uncertainty_label == "cv", sqrt(log(1 + obs$uncertainty^2)), obs$uncertainty) + expect_equal(res$se, se, label = name) + + p <- suppressMessages(plot_runs_test(dat, module = module, interactive = FALSE)) + expect_no_error(ggplot2::ggplot_build(p)) + } +}) + +test_that("index residuals keep every observation", { + two_seasons <- data.frame( + label = rep(c("index_observed", "index_predicted"), each = 6), + fleet = "Survey", + year = rep(rep(2001:2003, 2), 2), + season = rep(rep(1:2, each = 3), 2), + estimate = c(1, 2, 3, 10, 20, 30, 1.1, 2.1, 2.9, 11, 19, 31), + uncertainty = 0.2, + uncertainty_label = "se" + ) + res <- index_residuals(two_seasons, interactive = FALSE) + expect_setequal(unique(res$series), c("Survey (season 1)", "Survey (season 2)")) + # Exact copies are combined; values that cannot be told apart are an error + expect_equal(nrow(index_residuals(rbind(two_seasons, two_seasons), interactive = FALSE)), 6) + two_seasons$season <- NULL + expect_error(index_residuals(two_seasons, interactive = FALSE), "cannot be told apart") + + # Observations without a predicted value are reported + no_pred <- data.frame( + label = c(rep("index_observed", 3), rep("index_predicted", 2)), + fleet = "a", year = c(1:3, 1:2), estimate = c(1, 2, 3, 1, 2), + uncertainty = NA, uncertainty_label = NA + ) + expect_message(res <- index_residuals(no_pred, interactive = FALSE), "1 observation") + expect_equal(res$year, 1:2) + + expect_error( + index_residuals(list(stockplotr::example_data, stockplotr::example_data), interactive = FALSE), + "2 models" + ) + expect_error(index_residuals(stockplotr::example_data, fleet = "none", interactive = FALSE), "No index found") +}) + +test_that("plot_runs_test flags outliers and failed runs tests", { + p <- plot_runs_test(stockplotr::example_data, interactive = FALSE) + # The Triennial index in 2004 is outside the control limits + expect_equal(p$data$year[p$data$outside], 2004) + expect_equal(p$data$fleet[p$data$outside], "Triennial") + + p <- suppressMessages(plot_runs_test(stockplotr::pollock, interactive = FALSE)) + ctl <- p$layers[[1]]$data + expect_true(all(c("Passed", "Failed (p < 0.05)") %in% ctl$result)) + # Long fleet names are wrapped, not cut + expect_true(all(grepl("\n", ctl$panel))) +}) + +test_that("plot_runs_test saves a caption and alt text that describe the tests", { + rda <- fit_rda(plot_runs_test(stockplotr::example_data, interactive = FALSE, make_rda = TRUE), "runs_test") + expect_tags_filled(rda, "runs.test.") + expect_match(rda$alt_text, "passed for 2 of 2 indices") + expect_match(rda$alt_text, "One residual is outside the control limits: Triennial in 2004.") + # The results are recorded in key_quantities.csv + resid <- index_residuals(stockplotr::example_data, interactive = FALSE) + kqs <- rda$kqs[grepl("^runs[.]test[.]", names(rda$kqs))] + expect_setequal( + names(kqs), + paste0("runs.test.", c("start.year", "end.year", "n.indices", "n.passed", "failed", "p.values", "n.outside")) + ) + expect_equal(kqs[["runs.test.start.year"]], as.character(min(resid$year))) + expect_equal(kqs[["runs.test.n.indices"]], "2") + expect_equal(kqs[["runs.test.n.passed"]], "2") + expect_equal(kqs[["runs.test.failed"]], "none") + expect_equal(kqs[["runs.test.n.outside"]], "1") + expect_match(kqs[["runs.test.p.values"]], "^[^;]+: [0-9.]+; [^;]+: [0-9.]+$") + + # save_all_plots() saves to a figures folder other than the working directory + figs <- tempfile("figs") + on.exit(unlink(figs, recursive = TRUE), add = TRUE) + rda <- fit_rda( + plot_runs_test(stockplotr::example_data, interactive = FALSE, make_rda = TRUE, figures_dir = figs), + "runs_test", + dir = figs + ) + expect_tags_filled(rda, "runs.test.") +}) diff --git a/tests/testthat/test-save_all_plots.R b/tests/testthat/test-save_all_plots.R index 19eb6112..395960c8 100644 --- a/tests/testthat/test-save_all_plots.R +++ b/tests/testthat/test-save_all_plots.R @@ -10,6 +10,8 @@ fig_base_temp_files <- c( "catch_comp_figure.rda", "recruitment_deviations_figure.rda", "recruitment_figure.rda", + "joint_residuals_figure.rda", + "runs_test_figure.rda", "selectivity_figure.rda", "spawning_biomass_figure.rda", "spr_figure.rda", diff --git a/tests/testthat/test-utils_fit_captions.R b/tests/testthat/test-utils_fit_captions.R new file mode 100644 index 00000000..680882e0 --- /dev/null +++ b/tests/testthat/test-utils_fit_captions.R @@ -0,0 +1,50 @@ +test_that("caption text helpers format lists and panels", { + expect_equal(fit_list_text("A"), "A") + expect_equal(fit_list_text(c("A", "B")), "A and B") + expect_equal(fit_list_text(c("A", "B", "C", "A")), "A, B, and C") + expect_equal(fit_panel_text("a"), "A single panel for a.") + expect_equal(fit_panel_text(c("a", "b")), "A set of 2 panels, one for each fleet (a and b).") + expect_equal(fit_tag_value(c("costs $5", "a\\b")), "costs 5 ab") +}) + +test_that("fit_add_caption_rows adds missing rows to an existing captions file", { + dir <- tempfile("captions") + dir.create(dir) + old <- setwd(dir) + on.exit( + { + setwd(old) + unlink(dir, recursive = TRUE) + }, + add = TRUE + ) + template <- utils::read.csv( + system.file("resources", "captions_alt_text_template.csv", package = "stockplotr") + ) + utils::write.csv(template[template$label != "runs_test", ], "captions_alt_text.csv", row.names = FALSE) + fit_add_caption_rows("runs_test") + caps <- utils::read.csv("captions_alt_text.csv") + expect_equal(sum(caps$label == "runs_test"), 1) + expect_equal(nrow(caps), nrow(template)) + # Rows already present are not duplicated + fit_add_caption_rows("runs_test") + expect_equal(nrow(utils::read.csv("captions_alt_text.csv")), nrow(template)) +}) + +test_that("fit_add_kq_rows adds missing rows to an existing key quantities file", { + dir <- tempfile("kqs") + dir.create(dir) + on.exit(unlink(dir, recursive = TRUE), add = TRUE) + template <- utils::read.csv( + system.file("resources", "key_quantity_template.csv", package = "stockplotr") + ) + old <- template[!grepl("^runs[.]test[.]", template$key_quantity), ] + utils::write.csv(old, file.path(dir, "key_quantities.csv"), row.names = FALSE) + fit_add_kq_rows(c("runs.test.n.passed", "runs.test.failed"), dir = dir) + kqs <- utils::read.csv(file.path(dir, "key_quantities.csv")) + expect_equal(nrow(kqs), nrow(old) + 2) + expect_true(all(c("runs.test.n.passed", "runs.test.failed") %in% kqs$key_quantity)) + # Rows already present are not duplicated + fit_add_kq_rows("runs.test.n.passed", dir = dir) + expect_equal(nrow(utils::read.csv(file.path(dir, "key_quantities.csv"))), nrow(old) + 2) +})