diff --git a/.github/workflows/deploy-case-studies.yml b/.github/workflows/deploy-case-studies.yml new file mode 100644 index 0000000..cbdd920 --- /dev/null +++ b/.github/workflows/deploy-case-studies.yml @@ -0,0 +1,190 @@ +name: Deploy Case Studies + +on: + workflow_run: + workflows: ["Render Case Studies"] + types: + - completed + workflow_dispatch: + inputs: + run_id: + description: 'Run ID of Render Case Studies to get artifacts from (for testing)' + required: false + type: string + +permissions: + contents: write + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +jobs: + deploy: + if: > + github.event_name == 'workflow_dispatch' || github.event_name == 'push' || + ( + github.event.workflow_run.event == 'push' && + github.event.workflow_run.head_branch == 'main' && + (github.event.workflow_run.conclusion == 'success' || github.event.workflow_run.conclusion == 'failure') + ) + runs-on: ubuntu-latest + + steps: + - name: Check out repository + uses: actions/checkout@v4 + + - name: Download current run artifacts + uses: actions/download-artifact@v4 + with: + path: current-artifacts + github-token: ${{ secrets.GITHUB_TOKEN }} + run-id: ${{ inputs.run_id || github.event.workflow_run.id }} + + - name: Fetch existing gh-pages branch + shell: bash + run: | + git fetch origin gh-pages:gh-pages || true + + - name: Prepare site from previous deployment + shell: bash + run: | + rm -rf _site + mkdir -p _site + + if git show-ref --verify --quiet refs/heads/gh-pages; then + # The previous site is deployed to the root of the gh-pages branch. + # We extract it directly into the _site directory. + git archive gh-pages | tar -x -C _site/ + fi + + - name: Overlay successful renders + shell: bash + run: | + shopt -s nullglob + for artifact_dir in current-artifacts/*; do + if [ -d "$artifact_dir" ] && [[ "$(basename "$artifact_dir")" != status-* ]]; then + cp -R "$artifact_dir"/. _site/ + fi + done + + - name: Copy status files for README update + shell: bash + run: | + mkdir -p statuses + find current-artifacts -name '*.json' -path '*/status-*' -exec cp {} statuses/ \; + + - name: Setup R + uses: r-lib/actions/setup-r@v2 + + - name: Install jsonlite + run: Rscript -e 'install.packages("jsonlite", repos = "https://cloud.r-project.org")' + + - name: Update README table + shell: Rscript {0} + run: | + readme_file <- "README.md" + content_dir <- "content" + excluded <- c("setup.qmd", "acknowledgements.qmd", "advanced-features.qmd") + + parse_yaml <- function(filepath) { + lines <- readLines(filepath, warn = FALSE) + dash_lines <- which(lines == "---") + res <- list(stock = NA_character_, previous_model = "", features = "") + + if (length(dash_lines) >= 2) { + yaml_lines <- lines[(dash_lines[1] + 1):(dash_lines[2] - 1)] + for (line in yaml_lines) { + if (grepl("^stock:", line)) res$stock <- trimws(gsub("^stock:\\s*|\"|'", "", line)) + if (grepl("^previous_model:", line)) res$previous_model <- trimws(gsub("^previous_model:\\s*|\"|'", "", line)) + if (grepl("^features:", line)) res$features <- trimws(gsub("^features:\\s*|\"|'", "", line)) + } + } + res + } + + status_files <- list.files("statuses", pattern = "\\.json$", full.names = TRUE) + status_map <- list() + + if (length(status_files) > 0) { + for (sf in status_files) { + obj <- jsonlite::fromJSON(sf) + status_map[[obj$case_id]] <- obj$status + } + } + + files <- list.files(content_dir, pattern = "\\.qmd$", full.names = TRUE) + basenames <- basename(files) + valid_idx <- !(basenames %in% excluded) + files <- files[valid_idx] + basenames <- basenames[valid_idx] + + table_lines <- c( + "Stock | Previous Model | Status | Notable Features |", + "-- | -- | -- | --" + ) + + dir.create("badges", showWarnings = FALSE) + + for (i in seq_along(files)) { + file <- files[i] + case_id <- sub("\\.qmd$", "", basenames[i]) + meta <- parse_yaml(file) + + stock <- if (!is.na(meta$stock) && nzchar(meta$stock)) meta$stock else gsub("[-_]", " ", case_id) + status <- if (!is.null(status_map[[case_id]])) status_map[[case_id]] else "failing" + + badge_url <- if (status == "working") { + "https://img.shields.io/badge/Status-working-brightgreen" + } else { + "https://img.shields.io/badge/Status-failing-red" + } + + badge_path <- file.path("badges", paste0(case_id, ".svg")) + download.file(badge_url, destfile = badge_path, mode = "wb", quiet = TRUE) + + badge_md <- paste0("![Status](", badge_path, ")") + table_lines <- c(table_lines, paste(stock, meta$previous_model, badge_md, meta$features, sep = " | ")) + } + + readme <- readLines(readme_file, warn = FALSE) + start_idx <- which(trimws(readme) == "") + end_idx <- which(trimws(readme) == "") + + if (length(start_idx) > 0 && length(end_idx) > 0) { + new_readme <- c( + readme[1:start_idx], + table_lines, + readme[end_idx:length(readme)] + ) + writeLines(new_readme, readme_file) + } else { + stop("Could not find TABLE_START/TABLE_END markers in README.md") + } + + - name: Upload test artifacts + if: github.event_name == 'workflow_dispatch' + uses: actions/upload-artifact@v4 + with: + name: test-deploy-outputs + path: | + _site/ + README.md + + - name: Commit README update + if: github.event_name == 'workflow_run' + uses: peter-evans/create-pull-request@v6 + with: + commit-message: "chore: auto-update README case study table [skip ci]" + title: "🤖 Auto-update case study table" + body: "Automated update of the case study status table." + branch: "auto-update-readme-table" + base: "main" + + - name: Deploy merged site to gh-pages + if: github.event_name == 'workflow_run' + uses: peaceiris/actions-gh-pages@v4 + with: + github_token: ${{ secrets.GITHUB_TOKEN }} + publish_dir: ./_site + keep_files: true diff --git a/.github/workflows/render-case-studies.yml b/.github/workflows/render-case-studies.yml new file mode 100644 index 0000000..8297520 --- /dev/null +++ b/.github/workflows/render-case-studies.yml @@ -0,0 +1,204 @@ +name: Render Case Studies + +on: + push: + branches: [ main ] + paths: + - "content/**" + - ".github/workflows/render-case-studies.yml" + pull_request: + branches: [ main ] + paths: + - "content/**" + - ".github/workflows/render-case-studies.yml" + workflow_dispatch: + schedule: + - cron: "30 0 * * 0" + +permissions: + contents: read + actions: write + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +jobs: + discover: + runs-on: ubuntu-latest + outputs: + matrix: ${{ steps.set-matrix.outputs.matrix }} + steps: + - name: Check out repository + uses: actions/checkout@v4 + + - name: Setup R + uses: r-lib/actions/setup-r@v2 + + - name: Build matrix from content directory + id: set-matrix + shell: Rscript {0} + run: | + excluded <- c("setup.qmd", "acknowledgements.qmd", "advanced-features.qmd") + + files <- list.files("content", pattern = "\\.qmd$", full.names = FALSE) + files <- sort(files[!files %in% excluded]) + + include <- lapply(files, function(f) { + list( + qmd_path = file.path("content", f), + artifact_name = paste0(tools::file_path_sans_ext(f), "-site"), + case_id = tools::file_path_sans_ext(f) + ) + }) + + json_escape <- function(x) { + x <- gsub("\\\\", "\\\\\\\\", x) + x <- gsub('"', '\\"', x) + x + } + + entries <- vapply(include, function(x) { + sprintf( + '{"qmd_path":"%s","artifact_name":"%s","case_id":"%s"}', + json_escape(x$qmd_path), + json_escape(x$artifact_name), + json_escape(x$case_id) + ) + }, character(1)) + + matrix_json <- paste0('{"include":[', paste(entries, collapse = ","), ']}') + + cat("matrix=", matrix_json, "\n", sep = "", file = Sys.getenv("GITHUB_OUTPUT"), append = TRUE) + + render: + needs: discover + if: ${{ fromJson(needs.discover.outputs.matrix).include[0] != null }} + strategy: + fail-fast: false + matrix: ${{ fromJson(needs.discover.outputs.matrix) }} + runs-on: ubuntu-latest + + steps: + - name: Check out repository + uses: actions/checkout@v4 + + - name: Setup R + uses: r-lib/actions/setup-r@v2 + + - name: Install system dependencies + run: | + sudo apt-get update + sudo apt-get install -y \ + libcurl4-openssl-dev \ + libssl-dev \ + libxml2-dev \ + libfontconfig1-dev \ + libfreetype6-dev \ + libharfbuzz-dev \ + libfribidi-dev \ + libpng-dev \ + libtiff5-dev \ + libjpeg-dev \ + libcairo2-dev \ + libgit2-dev \ + libx11-dev \ + pandoc + + - name: Set up R package dependencies + uses: r-lib/actions/setup-r-dependencies@v2 + with: + pak-version: stable + cache: true + packages: | + cran::rmarkdown + cran::knitr + cran::jsonlite + cran::dplyr + cran::ggplot2 + cran::tidyr + cran::TMB + cran::remotes + cran::reshape2 + cran::Rcpp + cran::glue + cran::purrr + cran::ggridges + cran::yaml + extra-packages: | + github::nmfs-ost/stockplotr + github::stan-dev/cmdstanr + github::mjskay/tidybayes + github::stan-dev/shinystan + github::noaa-afsc/SparseNUTS + github::NOAA-FIMS/FIMS + github::r4ss/r4ss + + - name: Setup Quarto + uses: quarto-dev/quarto-actions/setup@v2 + + - name: Render case study + id: render + shell: bash + run: | + set -e + quarto render "${{ matrix.qmd_path }}" + + - name: Record status + if: always() + shell: bash + run: | + mkdir -p statuses + if [ "${{ steps.render.outcome }}" = "success" ]; then + status="working" + else + status="failing" + fi + + cat > "statuses/${{ matrix.case_id }}.json" <= 2) { - yaml_lines <- lines[(dash_lines[1] + 1):(dash_lines[2] - 1)] - for (line in yaml_lines) { - if (grepl("^stock:", line)) res$stock <- trimws(gsub("^stock:\\s*|\"|'", "", line)) - if (grepl("^previous_model:", line)) res$previous_model <- trimws(gsub("^previous_model:\\s*|\"|'", "", line)) - if (grepl("^features:", line)) res$features <- trimws(gsub("^features:\\s*|\"|'", "", line)) - } - } - return(res) - } - - # --- Initialize Table --- - table_lines <- c( - "Stock | Previous Model | Status | Notable Features |", - "-- | -- | -- | --" - ) - - # --- Process Quarto Files One at a Time --- - files <- list.files(content_dir, pattern = "\\.qmd$", full.names = TRUE) - basenames <- basename(files) - - # Filter out excluded files - valid_idx <- !(basenames %in% excluded) - files <- files[valid_idx] - basenames <- basenames[valid_idx] - - for (i in seq_along(files)) { - file <- files[i] - base <- gsub("\\.qmd$", "", basenames[i]) - - # Extract metadata - meta <- parse_yaml(file) - - # Fallback if stock isn't defined in YAML - stock <- ifelse(is.na(meta$stock), gsub("[-_]", " ", base), meta$stock) - - # Attempt to render the file - message(glue::glue("Testing {base}...")) - cmd <- glue::glue("quarto render {file}") - exit_code <- system(cmd) - - # Assign badge - if (exit_code == 0) { - message("✅ ", base, " rendered successfully.") - badge_url <- "https://img.shields.io/badge/Status-working-brightgreen" - } else { - message("❌ ", base, " failed to render.") - badge_url <- "https://img.shields.io/badge/Status-failing-red" - } - - # Download badge - badge_path <- file.path(badge_dir, paste0(base, ".svg")) - download.file(badge_url, destfile = badge_path, mode = "wb", quiet = TRUE) - - # Add row to table - badge_md <- glue::glue("![Status]({badge_dir}/{base}.svg)") - table_lines <- c(table_lines, glue::glue("{stock} | {meta$previous_model} | {badge_md} | {meta$features} |")) - } - - # --- Update README.md --- - readme <- readLines(readme_file, warn = FALSE) - start_idx <- which(trimws(readme) == "") - end_idx <- which(trimws(readme) == "") - - if (length(start_idx) > 0 && length(end_idx) > 0) { - new_readme <- c( - readme[1:start_idx], - table_lines, - readme[end_idx:length(readme)] - ) - writeLines(new_readme, readme_file) - message("README.md updated successfully!") - } else { - warning("Error: Could not find and markers in README.") - } - - - name: Deploy successfully rendered files (and not failed ones) to GitHub Pages - if: github.event_name == 'push' && github.ref == 'refs/heads/main' - uses: peaceiris/actions-gh-pages@v4 - with: - github_token: ${{ secrets.GITHUB_TOKEN }} - publish_dir: ./_site - keep_files: true - - - name: Create Pull Request - uses: peter-evans/create-pull-request@v8 - with: - # Added [skip ci] to prevent triggering workflows on merge - commit-message: "chore: auto-update case study table and badges [skip ci]" - title: "🤖 Auto-update case study table and badges" - body: "Automated updates to case study status badges and README table based on the latest Quarto renders." - branch: "auto-update-badges" - base: "main" diff --git a/.github/workflows/render.yml b/.github/workflows/render.yml index c751f12..539f4d8 100644 --- a/.github/workflows/render.yml +++ b/.github/workflows/render.yml @@ -1,64 +1,59 @@ +name: Render and Deploy Full Site +# for everything other than case studies + on: - pull_request: - branches: - - main - - dev push: - branches: - - dev + branches: [ main ] + paths: + - '!content/**' + - '!.github/workflows/render-case-studies.yml' + - 'README.md' workflow_dispatch: - schedule: - - cron: '30 0 * * 0' # runs At 00:00 on Sunday + schedule: + - cron: "0 2 * * 0" + +permissions: + contents: write + pages: write + id-token: write -name: Render +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true jobs: - build-deploy: + render: runs-on: ubuntu-latest - env: - GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} - GH_TOKEN: ${{ secrets.GITHUB_TOKEN }} steps: - name: Check out repository - uses: actions/checkout@v7 - - - name: Set up R (needed for Rmd) - id: setup-r - uses: r-lib/actions/setup-r@v2 + uses: actions/checkout@v4 - - name: Cache R package directories - uses: actions/cache@v6 - with: - path: | - ~/.cache/R - ~/R - key: ${{ runner.os }}-r-pak-${{ hashFiles('.github/workflows/render.yml') }} - restore-keys: | - ${{ runner.os }}-r-pak- + - name: Set up R + uses: r-lib/actions/setup-r@v2 - name: Install system dependencies run: | sudo apt-get update - sudo apt-get install -y \ - libcurl4-openssl-dev \ - libssl-dev \ - libxml2-dev \ - libv8-dev \ - libfontconfig1-dev \ - libfreetype6-dev \ - libharfbuzz-dev \ - libfribidi-dev \ - libpng-dev \ - libtiff5-dev \ - libjpeg-dev \ - libcairo2-dev \ - libgit2-dev \ - libnode-dev \ - libx11-dev \ - pandoc + sudo apt-get install -y \ + libcurl4-openssl-dev \ + libssl-dev \ + libxml2-dev \ + libv8-dev \ + libfontconfig1-dev \ + libfreetype6-dev \ + libharfbuzz-dev \ + libfribidi-dev \ + libpng-dev \ + libtiff5-dev \ + libjpeg-dev \ + libcairo2-dev \ + libgit2-dev \ + libnode-dev \ + libx11-dev \ + pandoc - - name: Set up R package dependencies (cached, via pak) + - name: Set up R package dependencies uses: r-lib/actions/setup-r-dependencies@v2 with: pak-version: stable @@ -88,15 +83,32 @@ jobs: - name: Set up Quarto uses: quarto-dev/quarto-actions/setup@v2 - env: - GH_TOKEN: ${{ secrets.GITHUB_TOKEN }} - - # add software dependencies here - - name: Render Quarto Project + - name: Render Quarto project uses: quarto-dev/quarto-actions/render@v2 - env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # this secret is always available for github actions - - name: Upload artifact - uses: actions/upload-pages-artifact@v5 + - name: Restore existing case-study pages from gh-pages + shell: bash + run: | + rm -rf /tmp/gh-pages + git clone --depth 1 --branch gh-pages "https://x-access-token:${{ secrets.GITHUB_TOKEN }}@github.com/${{ github.repository }}.git" /tmp/gh-pages || true + + if [ -d /tmp/gh-pages/content ]; then + mkdir -p _site/content + cp -R /tmp/gh-pages/content/. _site/content/ + fi + + - name: Upload pages artifact + uses: actions/upload-pages-artifact@v3 + with: + path: _site + + deploy: + needs: render + runs-on: ubuntu-latest + permissions: + pages: write + id-token: write + steps: + - name: Deploy to GitHub Pages + uses: actions/deploy-pages@v4 \ No newline at end of file diff --git a/content/AFSC-GOA-pollock.qmd b/content/AFSC-GOA-pollock.qmd index 421a28c..d7771ec 100644 --- a/content/AFSC-GOA-pollock.qmd +++ b/content/AFSC-GOA-pollock.qmd @@ -102,10 +102,10 @@ default_parameters <- FIMS::create_default_configurations( dplyr::rows_update( tibble::tibble( module_name = "Selectivity", - fleet_name = c("fleet1", "survey2", "survey6"), + fleet = c("fleet1", "survey2", "survey6"), module_type = "DoubleLogistic" ), - by = c("module_name", "fleet_name") + by = c("module_name", "fleet") ) |> FIMS::create_default_parameters( data = data_4_model @@ -114,7 +114,7 @@ default_parameters <- FIMS::create_default_configurations( dplyr::rows_update( tibble::tibble( module_name = "Selectivity", - fleet_name = "fleet1", + fleet = "fleet1", label = c( "inflection_point_asc", "slope_asc", "inflection_point_desc", "slope_desc" @@ -127,12 +127,12 @@ default_parameters <- FIMS::create_default_configurations( ), estimation_type = "fixed_effects" ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( module_name = "Selectivity", - fleet_name = "survey2", + fleet = "survey2", label = c( "inflection_point_asc", "slope_asc", "inflection_point_desc", "slope_desc" @@ -145,12 +145,12 @@ default_parameters <- FIMS::create_default_configurations( ), estimation_type = rep(c("fixed_effects", "constant"), each = 2) ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( module_name = "Selectivity", - fleet_name = "survey6", + fleet = "survey6", label = c( "inflection_point_asc", "slope_asc", "inflection_point_desc", "slope_desc" @@ -163,19 +163,19 @@ default_parameters <- FIMS::create_default_configurations( ), estimation_type = rep(c("constant", "fixed_effects"), each = 2) ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( module_name = "Selectivity", - fleet_name = "survey3", + fleet = "survey3", label = c("inflection_point", "slope"), value = c( parfinal$inf1_srv3, exp(parfinal$log_slp1_srv3) ) ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( @@ -238,22 +238,22 @@ default_parameters <- FIMS::create_default_configurations( dplyr::rows_update( tibble::tibble( module_name = "Fleet", - fleet_name = "fleet1", + fleet = "fleet1", time = years, label = "log_Fmort", value = log(pkfitfinal$rep$F) ), - by = c("module_name", "fleet_name", "label", "time") + by = c("module_name", "fleet", "label", "time") ) |> dplyr::rows_update( tibble::tibble( module_name = "Fleet", - fleet_name = c("survey2", "survey3", "survey6"), + fleet = c("survey2", "survey3", "survey6"), label = "log_q", value = c(parfinal$log_q2_mean, parfinal$log_q3_mean, parfinal$log_q6), estimation_type = "constant" ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) # Put it all together, creating the FIMS model and making the TMB fcn @@ -263,7 +263,7 @@ test_fit <- default_parameters |> FIMS::fit_fims(optimize = FALSE) FIMS::get_obj(test_fit)$report()$nll_components |> length() FIMS::get_data(data_4_model) |> - dplyr::filter(name == "fleet1", type == "age_comp") |> + dplyr::filter(fleet == "fleet1", type == "age_comp") |> print(n = 10) # Run the model with optimization fit <- default_parameters |> diff --git a/content/NEFSC-yellowtail.qmd b/content/NEFSC-yellowtail.qmd index 1d88cbb..242faa8 100644 --- a/content/NEFSC-yellowtail.qmd +++ b/content/NEFSC-yellowtail.qmd @@ -87,11 +87,11 @@ parameters <- create_default_configurations(data = data_4_model) |> # ASAP assumes Fmult devs = 0 dplyr::rows_update( tibble::tibble( - fleet_name = c("survey1", "survey2"), + fleet = c("survey1", "survey2"), label = "log_q", value = log(rdat$initial.guesses$q.year1.init) ), - by = c("fleet_name", "label") + by = c("fleet", "label") ) |> dplyr::rows_update( tibble::tibble( @@ -172,7 +172,7 @@ report <- get_report(fit) ```{r} #| label: comparison-plots -convert_output_asap <- function(report, fleet_names) { +convert_output_asap <- function(report, fleets) { years <- seq(report$parms$styr, report$parms$endyr) catch <- data.frame( module_name = "Fleet", @@ -248,7 +248,7 @@ convert_output_asap <- function(report, fleet_names) { dplyr::mutate( uncertainty = NA_real_, uncertainty_label = "se", - fleet = as.character(factor(module_id, labels = fleet_names)), + fleet = as.character(factor(module_id, labels = fleets)), era = "time" ) } @@ -316,7 +316,7 @@ print(comp_index) catch_results <- data.frame( observed = get_data(data_4_model) |> - dplyr::filter(type == "landings", name == "fishery") |> + dplyr::filter(type == "landings", fleet == "fishery") |> dplyr::pull(value), FIMS = report$landings_expected[[1]], ASAP = as.numeric(rdat$catch.pred[1,]) diff --git a/content/NWFSC-petrale.qmd b/content/NWFSC-petrale.qmd index 41c4a31..30c5621 100644 --- a/content/NWFSC-petrale.qmd +++ b/content/NWFSC-petrale.qmd @@ -15,35 +15,6 @@ format: ```{r} #| label: setup-objects common_name <- "Petrale sole" - - -# read SS3 input files from petrale sole assessment on GitHub -petrale_input <- r4ss::SS_read( - "https://raw.githubusercontent.com/pfmc-assessments/petrale/main/models/2023.a034.001/" -) -# add wtatage file to list of input files -petrale_input[["wtatage"]] <- r4ss::SS_readwtatage( - "https://raw.githubusercontent.com/pfmc-assessments/petrale/refs/heads/main/models/2023.a034.001/wtatage.ss_new" -) -# read SS3 output files from petrale sole assessment on GitHub -petrale_output <- r4ss::SS_output( - "https://raw.githubusercontent.com/pfmc-assessments/petrale/main/models/2023.a034.001", - printstats = FALSE, - verbose = FALSE, - covar = FALSE -) -output_petrale <- stockplotr::convert_output( - "https://raw.githubusercontent.com/pfmc-assessments/petrale/main/models/2023.a034.001/Report.sso", - model = "ss3" -) -# years <- seq(petrale_input$dat$styr, petrale_input$dat$endyr) -years <- seq(1896, petrale_input$dat$endyr) -n_years <- length(years) -# ages <- 0:petrale_input$dat$Nages # population ages in SS3, starts at age 0 -ages <- 1:17 # same as data bins -n_ages <- length(ages) -lengths <- petrale_output[["lbins"]] -nlengths <- length(lengths) ``` * R version: `r R_version` @@ -57,253 +28,49 @@ nlengths <- length(lengths) The model presented in this case study was changed substantially from the operational version and should not be considered reflective of the `r common_name` stock. These results are intended to be a demonstration and nothing more. -To get the operational model to more closely match a FIMS model the following changes were made: +To get the operational model to more closely match a FIMS model a simpler SS3 model was developed using code in https://github.com/NOAA-FIMS/case-studies/blob/petrale-with-new-r4ss-converter/content/R/NWFSC-petrale_simplify_assessment.R +These changes included * Remove data * Remove male lengths and ages * Remove discard fractions and discard comps * Simplify selectivity * Remove length-based retention functions - * Convert to age-based logistic (from length-based double-normal, fixed asymptotic) -* Remove priors - * Natural mortality - * Steepness -* Use female mean weight at age from the SS3 weight-at-age file + * Convert to age-based logistic (from asymptotic length-based double-normal with time blocks) +* Use time-invariant population-level female mean weight at age from the SS3 weight-at-age file for all fleets +* Use time-invariant maturity x fecundity at age from the SS3 weight-at-age file ## Setting up the data -```{r} -#| eval: false -#| label: local-production-model -#| echo: false - -## SIMPLIFY SS3 DATA FILE -dat <- petrale_input$dat - -# change dimensions -dat$Nsexes <- 1 # previously 2 -dat$Nages <- n_ages # previously 40 - -# mirroring selectivity and keeping both fishing fleets instead of aggregating -# catch aggregate landings across fleets -# dat$catch <- dat$catch |> -# dplyr::group_by(year) |> -# dplyr::summarize(catch = sum(catch), catch_se = mean(catch_se)) |> -# dplyr::mutate(seas = 1, fleet = 1, .after = year) |> -# as.data.frame() # write_fwf4() can't write tibles - -# filter indices (just include WCGBTS) -dat$CPUE <- dat$CPUE |> - dplyr::filter(index == 4) - -# remove discard info -dat$N_discard_fleets <- 0 -dat$discard_fleet_info <- NULL -dat$discard_data <- NULL - -# remove mean body weight -dat$use_meanbodywt <- 0 -dat$meanbodywt <- NULL -dat$DF_for_meanbodywt <- NULL - -# remove length comps -dat$use_lencomp <- 0 -dat$len_info <- NULL -dat$lencomp <- NULL - -# simplify ageing error matrix -dat$N_ageerror_definitions <- 1 -dat$ageerror <- dat$ageerror[1:2, paste0("age", 0:17)] - -# simplify age comps -dat$agecomp <- - dat$agecomp |> - # exclude fleets 2 and 3, and only use marginal data for fleet 4 - dplyr::filter(fleet %in% c(1, -4)) |> - # exclude male observations - dplyr::select(!dplyr::starts_with("m", ignore.case = FALSE)) |> - # assign observations to ageing error 1 - dplyr::mutate(Gender = 1, Ageerr = 1) |> - # remove redundant observations - # NOTE: this removes some age comp data because there - # were years with multiple observations from the same fleet - # due to multiple ageing error matrices - dplyr::distinct(year, fleet, .keep_all = TRUE) - -# Simplify SS3 control file -ctl <- petrale_input$ctl -ctl$Nages <- 17 # previously 40 -ctl$Nsexes <- 1 # previously 2 -# turn on empirical weight-at-age -# ctl$EmpiricalWAA <- 1 # previously zero -# remove male mortality and growth parameters -ctl$MG_parms <- ctl$MG_parms |> - dplyr::filter(!grepl("_Mal_", rownames(ctl$MG_parms))) -ctl$max_bias_adj <- -1 # set bias adjust = 1.0 for all years - -# filter catchability stuff to only fleet 4 -ctl$Q_options <- ctl$Q_options |> - dplyr::filter(fleet == 4) -ctl$Q_parms <- ctl$Q_parms |> - dplyr::filter(grepl("WCGBTS", rownames(ctl$Q_parms))) - -# add new age-based selex parameters -# age-based logistic (pattern 12) for fleets 1 and 4 -# have south fishery (fleet 2) mirror fleet 1 -# no parameters for fleet 3 (ignored index left in place to avoid renumbering fleets) -ctl$age_selex_types <- ctl$age_selex_types |> - dplyr::mutate( - Pattern = c(12, 15, 0, 12), - Special = c(0, 1, 0, 0) - ) # fleet 2 mirror fleet 1 - -# create 4 parameter rows -ctl$age_selex_parms <- cbind( - data.frame( - LO = 1, - HI = 10, - INIT = c(5, 2, 5, 2), - PRIOR = 0, - PR_SD = 99, - PR_type = 0, - PHASE = 2 - ), - matrix(0, nrow = 4, ncol = 7) -) -names(ctl$age_selex_parms)[8:14] <- names(ctl$size_selex_parms)[8:14] -rownames(ctl$age_selex_parms) <- paste0( - "AgeSel_", c("P_1_", "P_2_"), - c("fleet1", "fleet1", "fleet4", "fleet4") -) - -# remove length-based selectivity -ctl$size_selex_types <- ctl$size_selex_types |> - dplyr::mutate(Pattern = 0, Discard = 0, Male = 0) -ctl$size_selex_parms <- NULL -ctl$size_selex_parms_tv <- NULL - -# turn off variance adjustments -ctl$DoVar_adjust <- 0 -ctl$Variance_adjustment_list <- NULL -# use wtatage file -ctl$EmpiricalWAA <- 1 - -# Create wtatage.ss file -# wtatage <- SS_readwtatage( -# "c:/ss/Petrale/Petrale2023/petrale/models//2023.a034.001/wtatage.ss_new" -# ) -# wtatage for fleets -1:0 (population), and fleets 1:4 (real fleets) -wtatage <- matrix( - rep( - c( - 0.0010, # age 0 - 0.0148, 0.0617, 0.1449, 0.2570, 0.3876, 0.5260, 0.6640, 0.7957, 0.9175, - 1.0273, 1.1247, 1.2097, 1.2831, 1.3460, 1.3994, 1.4446, 1.4821 - ), - 6 - ), - nrow = 6, ncol = 18, byrow = TRUE -) -# wtatage for fleet = -2 (fecundity * maturity) -wtatage <- rbind( - c( - 0, 0, 0, 2.34337e-07, 6.7895e-06, 5.41665e-05, 0.000173396, 0.000338953, - 0.000514454, 0.000681685, 0.000834535, 0.000971472, 0.00109252, 0.00119837, - 0.00129006, 0.00136882, 0.00143602, 0.00149302 - ), - wtatage -) -names(wtatage) <- 0:17 -wtatage <- data.frame( - Yr = -1800, - Seas = 1, - Sex = 1, - Bio_Pattern = 1, - BirthSeas = 1, - Fleet = -2:4, - wtatage -) - -if (TRUE) { - # write modified data and control files - new_input <- petrale_input - new_input$dat <- dat - new_input$ctl <- ctl - new_input$wtatage <- wtatage - - # new_input$dir <- "c:/ss/Petrale/Petrale2023/petrale/models/2023.a050.001_FIMS_case-study" - new_input$dir <- "c:/ss/Petrale/Petrale2023/petrale/models/2023.a050.002_FIMS_case-study_wtatage" - # create directory if it doesn't exist - if (!dir.exists(new_input$dir)) { - dir.create(new_input$dir) - } - - r4ss::SS_write(new_input, dir = new_input$dir, overwrite = TRUE) - - r4ss::run( - new_input$dir, - show_in_console = TRUE, - skipfinished = FALSE - ) - p2 <- r4ss::SS_output(new_input$dir) -} -``` - ```{r} #| label: prepare-fims-data #| output: false #| warning: false - -# define the dimensions based on the range in the data -# for petrale, the north fleet catch begins in 1896, -# not the start year of the model 1876 (which is when the southern fleet catch begins, -# but for simplicity this petrale example only uses the northern catch) - -# fleet 4 used conditional age-at-length data with marginal observations -# entered as fleet == -4 (to exclude from likelihood due to redundancy) -# using only marginals for FIMS and exclude CAAL data by filtering out -# the fleet 4 age data -# similarly, some early ages were excluded from fleet 1 in the original model -# by assigning to fleet -1, these get filtered by get_ss3_data() - -# only include age comps with fleet = -4 or 1: -petrale_input$dat$agecomp <- petrale_input$dat$agecomp |> - dplyr::filter(fleet %in% c(-4, 1)) |> - dplyr::mutate(fleet = abs(fleet)) |> - dplyr::arrange(fleet, year, dplyr::desc(Nsamp)) |> - dplyr::distinct(fleet, year, .keep_all = TRUE) - -petrale_input$dat$lencomp <- petrale_input$dat$lencomp |> - dplyr::filter( - # filter out discard length-composition data - !(part == 1 & fleet == 1) - ) - -# convert SS3 data into FIMS format using function defined in the R directory -mydat <- get_ss3_data( - ss3_inputs = petrale_input, - fleets = c(1, 4), - ages = ages, - lengths = petrale_input[["dat"]][["lbin_vector"]] + +# convert SS3 data into FIMS format using function in r4ss and convert to FIMSFrame object +data_4_model <- r4ss::ss3_data_to_fims( + ss3_dir = "https://raw.githubusercontent.com/pfmc-assessments/petrale/main/models/2023.a050.003_FIMS_case-study_wtatage/", + ss_new = TRUE, + fleets = c(-4, 1, 2, 4) ) |> - # rename fleet4 as fleet2 (fleets 2 and 3 already removed above) - dplyr::mutate(name = dplyr::case_when( - name == "fleet1" ~ name, - name == "fleet4" ~ "survey1" # change fleet4 to survey1 to match yellowtail case-study - )) |> - dplyr::filter( - # filter out forecast years - timing %in% years - ) - - -# set up FIMS data objects -data_4_model <- FIMS::FIMSFrame(mydat |> - # Remove the length-composition data - # TODO: remove this as we create an age-to-length-conversion matrixt - dplyr::filter(type != "length_comp") |> - dplyr::select(-length) + FIMS::FIMSFrame() + + +# read SS3 output files from modified petrale sole assessment on GitHub +rep <- r4ss::SS_output( + "https://raw.githubusercontent.com/pfmc-assessments/petrale/main/models/2023.a050.003_FIMS_case-study_wtatage", + printstats = FALSE, + verbose = FALSE, + covar = FALSE ) + +# # commented out pending fix to https://github.com/nmfs-ost/stockplotr/pull/317 +# # alternative output processing using stockplotr +# output_petrale <- stockplotr::convert_output( +# "https://raw.githubusercontent.com/pfmc-assessments/petrale/main/models/2023.a050.003_FIMS_case-study_wtatage/Report.sso", +# model = "ss3" +# ) + ``` ## Run FIMS model @@ -312,220 +79,216 @@ data_4_model <- FIMS::FIMSFrame(mydat |> #| label: setup-model # lots of code below copied PIFSC-opakapaka.qmd -# rename SS3 output to avoid modifying lines copied from opakapaka model -rep <- petrale_output +years <- (data_4_model |> FIMS::get_start_year()):(data_4_model |> + FIMS::get_end_year()) +ages <- data_4_model |> FIMS::get_ages() # Create default parameters default_parameters <- FIMS::create_default_configurations( - data = data_4_model -) |> - FIMS::create_default_parameters( data = data_4_model - ) |> - tidyr::unnest(cols = data) - +) |> + FIMS::create_default_parameters( + data = data_4_model + ) |> + tidyr::unnest(cols = data) +# Note: the last value of the initial numbers at age is the first +# recruitment deviation +recdev_years <- default_parameters |> + dplyr::filter(module_name == "Recruitment" & label == "log_devs") |> + dplyr::pull(time) recdevs <- rep$recruitpars |> # $recruitpars is a subset of $parameters with Yr as an additional column - dplyr::filter(Yr %in% years[-1]) |> # filter out initial numbers parameters - dplyr::select(Yr, Value) - -init_naa <- rep$natage |> - dplyr::filter(Time == min(years), Sex == 1) |> # numbers at age in FIMS start year, not SS3 model - dplyr::select(paste(ages)) |> - as.numeric() # TODO: figure out if *1000 is needed -init_naa <- exp(rep$parameters |> dplyr::filter(grepl("R0", Label)) |> dplyr::pull(Init)) * 1000 * exp(-(ages - 1) * rep$parameters |> dplyr::filter(grepl("NatM.*Fem", Label)) |> dplyr::pull(Value)) -init_naa[n_ages] <- init_naa[n_ages] / rep$parameters |> dplyr::filter(grepl("NatM.*Fem", Label)) |> dplyr::pull(Value) # sum of infinite series + dplyr::filter(Yr %in% recdev_years) |> + dplyr::pull(Value) + +# initial numbers-at-age from SS3 output +init_naa_ss3 <- 1000 * + rep$natage |> + dplyr::filter(Time == min(years), Sex == 1) |> # numbers at age in FIMS start year, not SS3 model + dplyr::select(paste(ages)) |> + as.numeric() # TODO: figure out if *1000 is needed + +# alternative initial numbers-at-age based on exponential decay using female M +init_naa <- exp( + rep$parameters |> dplyr::filter(grepl("R0", Label)) |> dplyr::pull(Value) +) * + 1000 * + exp( + -(ages - 1) * + rep$parameters |> + dplyr::filter(grepl("NatM.*Fem", Label)) |> + dplyr::pull(Value) + ) +# init_naa calculation for plus group is based on sum of infinite series, +# which is equivalent to dividing the initial number at the last age by M +init_naa[length(ages)] <- init_naa[length(ages)] / + rep$parameters |> + dplyr::filter(grepl("NatM.*Fem", Label)) |> + dplyr::pull(Value) # sum of infinite series + +# ratio is close to 1, suggesting that approximation is reasonable +# range(init_naa / init_naa_ss3) +# # [1] 1.037574 1.228785 + +# function to get estimated parameter value from SS3 model +get_ss3_estimate <- function(string) { + rep$parameters |> + dplyr::filter(grepl( + pattern = string, + Label, + fixed = TRUE + )) |> + dplyr::pull(Value) +} parameters <- default_parameters |> - # SS3 model had length-based selectivity which leads to sex-specific - # age-based selectivity due to sexually-dimorphic growth. - # I didn't bother to calculate an age-based inflection point averaged over sexes - dplyr::rows_update( - tibble::tibble( - module_name = "Maturity", - label = c("inflection_point", "slope"), - value = c(20.5, 1.8), # arbitrary guess of slope - ), - by = c("module_name", "label") - ) |> - dplyr::rows_update( - tibble::tibble( - module_name = "Selectivity", - fleet_name = "fleet1", - label = c("inflection_point", "slope"), - value = c(10, 2), - estimation_type = "constant" - ), - by = c("module_name", "fleet_name", "label") - ) |> - dplyr::rows_update( - tibble::tibble( - module_name = "Selectivity", - fleet_name = "survey1", - label = c("inflection_point", "slope"), - value = c(6, 2), - estimation_type = "constant" - ), - by = c("module_name", "fleet_name", "label") - ) |> - dplyr::rows_update( - tibble::tibble( - module_name = "Recruitment", - label = c("log_rzero", "logit_steep", "log_sd"), - value = c( - petrale_input$ctl$SR_parms["SR_LN(R0)", "INIT"], - FIMS::logit(0.2, 1.0, petrale_input$ctl$SR_parms["SR_BH_steep", "INIT"]), - 0.5 - ) - ), - by = c("module_name", "label") - ) |> - dplyr::rows_update( - tibble::tibble( - module_name = "Recruitment", - label = "log_devs", - time = (FIMS::get_start_year(data_4_model) + 1):FIMS::get_end_year(data_4_model), - # The last value of the initial numbers at age is the first - # recruitment deviation - value = recdevs$Value, - ), - by = c("module_name", "label", "time") - ) |> - dplyr::rows_update( - tibble::tibble( - module_name = "Population", - label = "log_init_naa", - age = FIMS::get_ages(data_4_model), - value = log(init_naa), - estimation_type = "fixed_effects" - ), - by = c("module_name", "label", "age") - ) |> - dplyr::rows_update( - tibble::tibble( - module_name = "Population", - label = "log_M", - value = log(rep$parameters["NatM_uniform_Fem_GP_1", "Value"]), - estimation_type = "constant" - ), - by = c("module_name", "label") - ) |> - dplyr::rows_update( - tibble::tibble( - module_name = "Fleet", - fleet_name = "fleet1", - time = years, - label = "log_Fmort", - estimation_type = "constant", - value = rep[["exploitation"]] |> - dplyr::filter(Yr %in% years) |> - dplyr::mutate(North = ifelse(North == 0, -999, log(North))) |> - dplyr::pull(North) - ), - by = c("module_name", "fleet_name", "label", "time") - ) + # dplyr::rows_update( + # tibble::tibble( + # module_name = "Maturity", + # label = c("inflection_point", "slope"), + # value = c(7.0, 1.0), # maturity based on the derived age-based pattern from original length-based SS3 model, arbitrary guess of slope + # ), + # by = c("module_name", "label") + # ) |> + # dplyr::rows_update( + # tibble::tibble( + # module_name = "Selectivity", + # fleet = "North", + # label = c("inflection_point", "slope"), + # value = c(10, 2), + # estimation_type = "constant" + # ), + # by = c("module_name", "fleet", "label") + # ) |> + # dplyr::rows_update( + # tibble::tibble( + # module_name = "Selectivity", + # fleet = "South", + # label = c("inflection_point", "slope"), + # value = c(10, 2), + # estimation_type = "constant" + # ), + # by = c("module_name", "fleet", "label") + # ) |> + # dplyr::rows_update( + # tibble::tibble( + # module_name = "Selectivity", + # fleet = "WCGBTS", + # label = c("inflection_point", "slope"), + # value = c(6, 2), + # estimation_type = "constant" + # ), + # by = c("module_name", "fleet", "label") + # ) |> + + # set initial values for stock-recruitment parameters based on SS3 output + dplyr::rows_update( + tibble::tibble( + module_name = "Recruitment", + label = c("log_rzero", "logit_steep", "log_sd"), + value = c( + log(1000) + # adding log(1000) to account for R0 units in 1000s in SS3 + get_ss3_estimate("SR_LN(R0)"), + FIMS::logit( + 0.2, + 1.0, + get_ss3_estimate("SR_BH_steep") + ), + get_ss3_estimate("SR_sigmaR") + ) + ), + by = c("module_name", "label") + ) |> + # set initial values for recruitment deviations based on SS3 output, + # starting with the second year of the SS3 model since the first recruitment deviation + # corresponds to the initial numbers at age + dplyr::rows_update( + tibble::tibble( + module_name = "Recruitment", + label = "log_devs", + time = recdev_years, + value = recdevs, + ), + by = c("module_name", "label", "time") + ) |> + # set initial value for log_q based on SS3 output + dplyr::rows_update( + tibble::tibble( + fleet = "WCGBTS", + label = "log_q", + value = get_ss3_estimate("LnQ_base_WCGBTS(4)") + ), + by = c("fleet", "label") + ) |> + # set initial numbers at age based on SS3 output calculated above + dplyr::rows_update( + tibble::tibble( + module_name = "Population", + label = "log_init_naa", + age = FIMS::get_ages(data_4_model), + value = log(init_naa), + estimation_type = "fixed_effects" + ), + by = c("module_name", "label", "age") + ) |> + # set initial M based on SS3 output + dplyr::rows_update( + tibble::tibble( + module_name = "Population", + label = "log_M", + value = log(get_ss3_estimate("NatM_uniform_Fem_GP_1")), + estimation_type = "constant" + ), + by = c("module_name", "label") + ) + +# confirm that all columns except "value" and "estimation_type" are identical between parameters and default_parameters +if (!identical( + parameters |> + dplyr::select(-value, -estimation_type), + default_parameters |> + dplyr::select(-value, -estimation_type) +)) { + cli::cli_alert_danger("Columns other than 'value' and 'estimation_type' are not identical between parameters and default_parameters") +} +# figure out if any rows have different estimate_type +cli::cli_alert_info("Rows with different estimation_type between default_parameters parameters:") +diff_estimation_type <- parameters |> + dplyr::filter(estimation_type != default_parameters$estimation_type) |> + dplyr::select(module_name, fleet, label, estimation_type) + +if (nrow(diff_estimation_type) == 0) { + cli::cli_alert_success("No rows with different estimation_type.") +} else { + print(diff_estimation_type, n = Inf, width = Inf) +} +# parameters which have been updated from defaults +cli::cli_alert_info("Rows with updated initial values from default_parameters to parameters:") +updated_parameter_counts <- parameters |> + dplyr::filter(value != default_parameters$value) |> + dplyr::select(module_name, fleet, label, value) |> + dplyr::count(module_name, label) + +if (nrow(updated_parameter_counts) == 0) { + cli::cli_alert_success("No rows with updated initial values.") +} else { + print(updated_parameter_counts, n = Inf, width = Inf) +} # Run the model without optimization to help ensure a viable model test_fit <- parameters |> - FIMS::initialize_fims(data = data_4_model) |> - FIMS::fit_fims(optimize = FALSE) + FIMS::initialize_fims(data = data_4_model) |> + FIMS::fit_fims(optimize = FALSE) # Run the model with optimization fit <- parameters |> - FIMS::initialize_fims(data = data_4_model) |> - FIMS::fit_fims(optimize = TRUE, get_sd = FALSE) + FIMS::initialize_fims(data = data_4_model) |> + FIMS::fit_fims(optimize = TRUE, get_sd = FALSE) +fit_est <- fit # backup in case we want to set fit <- test_fit ``` ## Plotting results -```{r} -#| label: comparison-plots -# gather index fit info -index_results <- data.frame( - observed = FIMS::m_index(data_4_model, "survey1"), - expected = FIMS::get_report(fit)[["index_expected"]][[2]] -) |> - dplyr::mutate(year = years) |> - dplyr::filter(observed > 0) # filter out -999 rows - -# plot index fit -ggplot2::ggplot(index_results, ggplot2::aes(x = year, y = observed)) + - ggplot2::geom_point() + - ggplot2::xlab("Year") + - ggplot2::ylab("Index (mt)") + - ggplot2::geom_line(ggplot2::aes(x = year, y = expected), color = "blue") + - ggplot2::theme_bw() - - -# gather catch info -catch_results <- data.frame( - year = years, - observed = FIMS::m_landings(data_4_model, fleet = "fleet1"), - expected = FIMS::get_report(fit)[["landings_expected"]][[1]] -) - -# plot catch fit -ggplot2::ggplot(catch_results, ggplot2::aes(x = year, y = observed)) + - ggplot2::geom_point() + - #ggplot2::aes(color = fleet)) + - ggplot2::xlab("Year") + - ggplot2::ylab("Catch (mt)") + - ggplot2::geom_line(ggplot2::aes(x = year, y = expected)) + #, color = fleet)) + - ggplot2::theme_bw() - -# gather biomass info from FIMS and SS3 -biomass <- rep$timeseries |> - dplyr::select(Yr, SpawnBio, Bio_all) |> - dplyr::filter(Yr %in% years) |> - dplyr::rename( - "SS3_SpawnBio" = "SpawnBio", - "SS3_Bio" = "Bio_all" - ) |> - dplyr::mutate( - FIMS_SpawnBio = FIMS::get_report(fit)[["spawning_biomass"]][[1]][-1] , - FIMS_Bio = FIMS::get_report(fit)[["biomass"]][[1]][-1] - ) |> ##CHECK: Is FIMS ssb reporting n_years+1 or initial year-1? - tidyr::pivot_longer(cols = -Yr) |> - tidyr::separate_wider_delim( - cols = "name", - delim = "_", names = c("Model", "Type") - ) - -# plot comparison of spawning biomass and total biomass time series -ggplot2::ggplot(biomass, ggplot2::aes(x = Yr, y = value)) + - ggplot2::geom_line(ggplot2::aes(color = Model)) + - ggplot2::xlab("Year") + - ggplot2::ylab("") + - ggplot2::facet_wrap(~Type, scales = "free_y") + - ggplot2::theme_bw() - -# gather info on recruitment -recruits <- rep$recruit |> - dplyr::filter(Yr %in% years) |> - dplyr::select(Yr, exp_recr, raw_dev) |> - dplyr::rename("SS3_recruit" = "exp_recr", - "SS3_recdev" = "raw_dev", - "Year" = "Yr") |> - dplyr::mutate(FIMS_recruit = c(FIMS::get_report(fit)[["expected_recruitment"]][[1]][-1]), - FIMS_recdev = c( - NA, - FIMS::get_estimates(fit) |> - dplyr::filter( - module_name == "Recruitment", - label == "log_devs" - ) |> - dplyr::pull(estimated) - ), - SS3_recruit = SS3_recruit * 1000) |> - tidyr::pivot_longer(cols = -Year) |> - tidyr::separate_wider_delim(cols = "name", delim = "_", names = c("Model", "Type")) - -# plot recruit time series -ggplot2::ggplot(recruits, ggplot2::aes(x = Year, y = value, color = Model)) + - ggplot2::geom_line() + - ggplot2::facet_wrap(~Type, scales = "free_y") + - ggplot2::theme_bw() - -``` +TODO: Add plotting code for the FIMS model results ```{r} #| label: cleanup diff --git a/content/PIFSC-opakapaka.qmd b/content/PIFSC-opakapaka.qmd index 9306443..7f1ecfe 100644 --- a/content/PIFSC-opakapaka.qmd +++ b/content/PIFSC-opakapaka.qmd @@ -1,8 +1,5 @@ --- -title: PIFSC Opakapaka Case Study -stock: "PIFSC Opakapaka" -previous_model: "SS3" -features: "Age-to-length conversion matrix" +title: PIFSC ʻŌpakapaka Case Study format: html: code-fold: true @@ -40,323 +37,302 @@ To get the operational model to more closely match a FIMS model the following ch #| eval: false #| label: local-production-model #| echo: false - # read in data and control files of original model # Locally this is # C:/Users/Megumi.Oshima/Documents/Opaka-FIMS-Case-Study/Model/01_original_model -opaka_mod_dir <- file.path( - getwd(), - "..", - "Opaka-FIMS-Case-Study", - "Model", - "01_original_model" -) -opaka_dat <- r4ss::SS_readdat_3.30(file.path(opaka_mod_dir, "data.ss")) -opaka_ctl <- r4ss::SS_readctl_3.30( - file.path(opaka_mod_dir, "control.ss"), - datlist = file.path(opaka_mod_dir, "data.ss") -) - -# create directory for new simplified model -opaka_length_dir <- file.path(opaka_mod_dir, "..", "09_case_study_lengths") -dir.create(opaka_length_dir) - -# remove size freq data -opaka_dat$N_sizefreq_methods_rd <- 0 -opaka_dat$N_sizefreq_methods <- NULL -opaka_dat$nbins_per_method <- NULL -opaka_dat$units_per_method <- NULL -opaka_dat$scale_per_method <- NULL -opaka_dat$mincomp_per_method <- NULL -opaka_dat$Nobs_per_method <- NULL -opaka_dat$Comp_Error_per_method <- NULL -opaka_dat$ParmSelect_per_method <- NULL -opaka_dat$sizefreq_bins_list <- NULL -opaka_dat$sizefreq_data_list <- NULL - -# remove super periods for length comp data -len_dat_original <- read.csv( - file.path(opaka_mod_dir, "..", "..", "Data", "Opaka_len_data.csv") -) -opaka_dat$lencomp <- len_dat_original - -# remove dirichlet weighting for length comps -opaka_dat$len_info$CompError <- 0 -opaka_dat$len_info$ParmSelect <- 0 - -# remove initial F estimation -opaka_dat$catch[1,'catch'] <- 0 -# add agecomp dummy data -opaka_dat$N_agebins <- 21 -opaka_dat$agebin_vector <- seq(1, 21) -opaka_dat$N_ageerror_definitions <- 1 -opaka_dat$ageerror <- rbind( + +# code to modify original model and create bootstrap data +# this creates the .RDS file which is included with the case studies but is not run when rendering the case study. The code is included here for transparency and reproducibility but is not intended to be run by users of the case study. +if (FALSE) { + # local directory to write modified data and control files to and run SS3 + local_dir <- tempdir() + + opaka_inputs_original <- r4ss::SS_read( + "https://raw.githubusercontent.com/MOshima-PIFSC/Opaka-FIMS-Case-Study/refs/heads/main/Model/01_original_model/", + ss_new = TRUE + ) + + opaka_dat <- opaka_inputs_original$dat + opaka_ctl <- opaka_inputs_original$ctl + + # create directory for new simplified model + opaka_length_dir <- file.path(local_dir, "09_case_study_lengths") + dir.create(opaka_length_dir) + + # changes to data file + + # remove size freq data + opaka_dat$N_sizefreq_methods_rd <- 0 + opaka_dat$N_sizefreq_methods <- NULL + opaka_dat$nbins_per_method <- NULL + opaka_dat$units_per_method <- NULL + opaka_dat$scale_per_method <- NULL + opaka_dat$mincomp_per_method <- NULL + opaka_dat$Nobs_per_method <- NULL + opaka_dat$Comp_Error_per_method <- NULL + opaka_dat$ParmSelect_per_method <- NULL + opaka_dat$sizefreq_bins_list <- NULL + opaka_dat$sizefreq_data_list <- NULL + + # remove super periods for length comp data + len_dat_original <- read.csv( + "https://raw.githubusercontent.com/MOshima-PIFSC/Opaka-FIMS-Case-Study/refs/heads/main/Data/Opaka_len_data.csv" + ) + opaka_dat$lencomp <- len_dat_original + + # remove dirichlet weighting for length comps + opaka_dat$len_info$CompError <- 0 + opaka_dat$len_info$ParmSelect <- 0 + + # remove initial F estimation + opaka_dat$catch[1, "catch"] <- 0 + # add agecomp dummy data + opaka_dat$N_agebins <- 21 + opaka_dat$agebin_vector <- seq(1, 21) + opaka_dat$N_ageerror_definitions <- 1 + opaka_dat$ageerror <- rbind( seq(0.5, 43.5), rep(.01, 44) -) -opaka_dat$age_info <- data.frame( - mintailcomp <- rep(0, 4), - addtocomp = 1e-7, - combine_M_F = 1, - CompressBins = 0, - CompError = 0, - ParmSelect = 0, - minsamplesize = 1 -) -opaka_dat$Lbin_method <- 1 -agecomp_info <- data.frame( - year = rep(seq(2017, 2023), 2), - month = 1, - fleet = rep(c(2, 4), each = 7), - sex = 0, - part = 0, - ageerr = 1, - Lbin_lo = -1, - Lbin_hi = -1, - Nsamp = 10 -) -dummy_agecomp <- as.data.frame(matrix( - data = 1, - nrow = nrow(agecomp_info), - ncol = length(opaka_dat$agebin_vector) -)) -colnames(dummy_agecomp) <- paste0("a", opaka_dat$agebin_vector) -opaka_dat$agecomp <- cbind(agecomp_info, dummy_agecomp) - -r4ss::SS_writedat_3.30( - datlist = opaka_dat, - outfile = file.path(opaka_length_dir, "data.ss"), - overwrite = TRUE -) - -# remove growth platoon -opaka_ctl$N_platoon <- 1 -opaka_ctl$sd_ratio <- NULL -opaka_ctl$submorphdist <- NULL - -# remove intial F estimation -opaka_ctl$init_F <- NULL - -# remove extra SE parameter -opaka_ctl$Q_options[1, "extra_se"] <- 0 -opaka_ctl$Q_parms <- opaka_ctl$Q_parms[-2, ] -opaka_ctl$Variance_adjustment_list <- NULL -opaka_ctl$DoVar_adjust <- 0 -opaka_ctl$sd_offset <- 0 - -# remove dirichlet weighting parameter lines -opaka_ctl$dirichlet_parms <- NULL - -# fix commercial selectivity -opaka_ctl$size_selex_parms$PHASE[1:2] <- -2 - -# add age selectivity -opaka_ctl$age_selex_types <- data.frame( - Pattern = rep(12, 4), - Discard = 0, - Male = 0, - Special = 0 -) - -# control file wouldn't write when age_selex_params are manually added -# opaka_ctl$age_selex_parms <- data.frame( -# "LO" = c(0, -5, 0, 0, 0, -10, 0, -20), -# "HI" = c(40, 50, 40, 40, 60, 60, 10, 50), -# "INIT" = c(1.81975, 0.0093046, 1, 3, 1.97182, 0.00040, 1.29111, 0.00115), -# "PRIOR" = c(5, 6, 5, 6, 5, 6, 2, .5), -# "PR_SD" = c(99, 99, 99, 99, 99, 99, 5, 2), -# "PR_type" = 0, -# "PHASE" = c(-2, -2, -2, -2, -99, -99, -2, -2), -# "env-var" = 0, -# "use_dev" = 0, -# "dev_mnyr" = 0, -# "dev_mxyr" = 0, -# "dev_PH" = 0, -# "Block" = 0, -# "Block_Fxn" = 0 -# ) - -age_ctl <- r4ss::SS_readctl_3.30( - file.path(opaka_mod_dir, "..", "03_age_comps", "control.ss_new"), - datlist = file.path(opaka_mod_dir, "..", "03_age_comps", "data.ss") -) -age_selex_params <- age_ctl$age_selex_parms -opaka_ctl$age_selex_parms <- age_selex_params -opaka_ctl$age_selex_parms$PHASE <- -2 - -r4ss::SS_writectl_3.30( - opaka_ctl, - outfile = file.path(opaka_length_dir, "control.ss"), - overwrite = TRUE -) - -ss_files <- c("forecast.ss", "starter.ss", "ss_opt_win.exe") -file.copy(file.path(opaka_mod_dir, ss_files), opaka_length_dir) - -# create a bootstrap data file to get age comp data -start <- r4ss::SS_readstarter(file.path(opaka_length_dir, "starter.ss")) -start$N_bootstraps <- 3 -r4ss::SS_writestarter(start, dir = opaka_length_dir, overwrite = TRUE) - -# run SS3 -r4ss::run(dir = opaka_length_dir, exe = "ss_opt_win.exe", skipfinished = FALSE) - -file.copy( - file.path(opaka_length_dir, "data_boot_001.ss"), - file.path(opaka_length_dir, "data.ss"), - overwrite = TRUE -) -start <- r4ss::SS_readstarter(file.path(opaka_length_dir, "starter.ss")) -start$N_bootstraps <- 1 -r4ss::SS_writestarter(start, dir = opaka_length_dir, overwrite = T) - -# run SS3 -r4ss::run(dir = opaka_length_dir, exe = "ss_opt_win.exe", skipfinished = F) - -# check model -rep <- r4ss::SS_output(dir = opaka_length_dir) -SS_plots(rep) - -# package up data, control and rep file for using in FIMS -rm("opaka_dat") -rm("opaka_ctl") -opaka_dat <- r4ss::SS_readdat_3.30(file.path(opaka_length_dir, "data.ss")) -opaka_ctl <- r4ss::SS_readctl_3.30( - file.path(opaka_length_dir, "control.ss"), - datlist = file.path(opaka_length_dir, "data.ss") -) -save( - list = c("opaka_dat", "opaka_ctl", "rep"), - file = file.path(opaka_length_dir, "opaka_length.RDS") -) + ) + opaka_dat$age_info <- data.frame( + mintailcomp = rep(0, 4), + addtocomp = 1e-7, + combine_M_F = 1, + CompressBins = 0, + CompError = 0, + ParmSelect = 0, + minsamplesize = 1 + ) + opaka_dat$Lbin_method <- 1 + agecomp_info <- data.frame( + year = rep(seq(2017, 2023), 2), + month = 1, + fleet = rep(c(2, 4), each = 7), + sex = 0, + part = 0, + ageerr = 1, + Lbin_lo = -1, + Lbin_hi = -1, + Nsamp = 10 + ) + dummy_agecomp <- as.data.frame(matrix( + data = 1, + nrow = nrow(agecomp_info), + ncol = length(opaka_dat$agebin_vector) + )) + colnames(dummy_agecomp) <- paste0("a", opaka_dat$agebin_vector) + opaka_dat$agecomp <- cbind(agecomp_info, dummy_agecomp) + + # changes to control file + + # remove growth platoon + opaka_ctl$N_platoon <- 1 + opaka_ctl$sd_ratio <- NULL + opaka_ctl$submorphdist <- NULL + + # remove initial F estimation + opaka_ctl$init_F <- NULL + + # remove extra SE parameter + opaka_ctl$Q_options[1, "extra_se"] <- 0 + opaka_ctl$Q_parms <- opaka_ctl$Q_parms[-2, ] + opaka_ctl$Variance_adjustment_list <- NULL + opaka_ctl$DoVar_adjust <- 0 + opaka_ctl$sd_offset <- 0 + + # remove dirichlet weighting parameter lines + opaka_ctl$dirichlet_parms <- NULL + + # fix commercial selectivity + opaka_ctl$size_selex_parms$PHASE[1:2] <- -2 + + # add age selectivity + opaka_ctl$age_selex_types <- data.frame( + Pattern = rep(12, 4), + Discard = 0, + Male = 0, + Special = 0 + ) + + # control file wouldn't write when age_selex_params are manually added + # opaka_ctl$age_selex_parms <- data.frame( + # "LO" = c(0, -5, 0, 0, 0, -10, 0, -20), + # "HI" = c(40, 50, 40, 40, 60, 60, 10, 50), + # "INIT" = c(1.81975, 0.0093046, 1, 3, 1.97182, 0.00040, 1.29111, 0.00115), + # "PRIOR" = c(5, 6, 5, 6, 5, 6, 2, .5), + # "PR_SD" = c(99, 99, 99, 99, 99, 99, 5, 2), + # "PR_type" = 0, + # "PHASE" = c(-2, -2, -2, -2, -99, -99, -2, -2), + # "env-var" = 0, + # "use_dev" = 0, + # "dev_mnyr" = 0, + # "dev_mxyr" = 0, + # "dev_PH" = 0, + # "Block" = 0, + # "Block_Fxn" = 0 + # ) + + age_ctl <- r4ss::SS_read( + "https://raw.githubusercontent.com/MOshima-PIFSC/Opaka-FIMS-Case-Study/refs/heads/main/Model/03_age_comps/", + ss_new = TRUE + )$ctl + + age_selex_params <- age_ctl$age_selex_parms + opaka_ctl$age_selex_parms <- age_selex_params + opaka_ctl$age_selex_parms$PHASE <- -2 + + # gather modified elements into new list of inputs + opaka_inputs_modified <- opaka_inputs_original + opaka_inputs_modified$dat <- opaka_dat + opaka_inputs_modified$ctl <- opaka_ctl + # create a bootstrap data file to get age comp data + opaka_inputs_modified$start$N_bootstraps <- 3 + + # write modified input files to local directory to run SS3 + r4ss::SS_write( + opaka_inputs_modified, + dir = opaka_length_dir, + overwrite = TRUE + ) + + + ## specify SS3 exe name and get exe from repo: + exe <- "ss3" + r4ss::get_ss3_exe(dir = opaka_length_dir, version = "v3.30.21") + ## run SS3 + r4ss::run(dir = opaka_length_dir, exe = exe, skipfinished = FALSE) + + file.copy( + file.path(opaka_length_dir, "data_boot_001.ss"), + file.path(opaka_length_dir, "data.ss"), + overwrite = TRUE + ) + + # run SS3 using bootstrap data + r4ss::run(dir = opaka_length_dir, exe = exe, skipfinished = FALSE) + + # check model + rep <- r4ss::SS_output(dir = opaka_length_dir) + r4ss::SS_plots(rep) + + # package up data, control and rep file for using in FIMS + rm("opaka_dat") + rm("opaka_ctl") + + # read input files again (to get the bootstrap data and the wtatage.ss_new) + opaka_inputs_bootstrap <- r4ss::SS_read(opaka_length_dir, read_wtatage = TRUE) + # save to model directory + save( + list = c("opaka_inputs_bootstrap", "rep"), + file = file.path(opaka_length_dir, "opaka_model.RDS") + ) + # copy to data directory for case studies repository + file.copy( + file.path(opaka_length_dir, "opaka_model.RDS"), + file.path(data_directory, "opaka_model.RDS"), + overwrite = TRUE + ) +} # end if (FALSE) for code used to prepare .RDS file ``` ```{r} #| label: prepare-fims-data #| output: false #| warning: false +# load opaka_inputs_bootstrap and rep objects with model input and output load(file.path(data_directory, "opaka_model.RDS")) include_age_comps <- FALSE -years <- seq(opaka_dat$styr, opaka_dat$endyr) +years <- seq(opaka_inputs_bootstrap$dat$styr, opaka_inputs_bootstrap$dat$endyr) +alk_years <- c(years, max(years) + 1) n_years <- length(years) # the number of years which we have data for. ages <- seq(1, 21) # age vector. n_ages <- length(ages) # the number of age groups. -comp_lengths <- opaka_dat$lbin_vector # length vector. +comp_lengths <- opaka_inputs_bootstrap$dat$lbin_vector # length vector. nlengths <- length(comp_lengths) # the number of length bins. -opaka_dat_fims <- get_ss3_data( - list( - dat = opaka_dat, - ctl = opaka_ctl, - start = list(), - fore = list(), - wtatage = rep[["wtatage"]] - ), - fleets = c(1,2,3), - ages = ages, +opaka_dat_fims_raw <- r4ss::ss3_data_to_fims( + ss3_inputs = opaka_inputs_bootstrap, + ss3_output = rep, + fleets = c(1, 2, 3), + maxage = max(ages), lengths = comp_lengths ) |> dplyr::filter(type != "age_comp") -## age to length conversion matrix -# Growth function values to create age to length conversion matrix from model -#comparison project -mg_pars <- rep$parameters |> -dplyr::filter(stringr::str_detect(Label, "_GP_")) -Linf <- mg_pars$Value[3] -K <- mg_pars$Value[4] -a0 <- -0.29 -amax <- 21 -cv <- mg_pars$Value[5] - -L2Wa <- mg_pars$Value[7] -L2Wb <- mg_pars$Value[8] - -AtoL <- function(a,Linf,K,a_0){ - L <- Linf*(1-exp(-K*(a-a_0))) - } - -ages <- 1:amax -len_bins <- comp_lengths - -#Create length at age conversion matrix and fill proportions using above -#growth parameters -length_age_conversion <- matrix(NA,nrow=length(ages),ncol=length(len_bins)) -for(i in seq_along(ages)){ - #Calculate mean length at age to spread lengths around - mean_length <- AtoL(ages[i],Linf,K,a0) - #mean_length <- AtoLSchnute(ages[i],L1,L2,a1,a2,Ks) - #Calculate the cumulative proportion shorter than each composition length - temp_len_probs<-pnorm(q=len_bins,mean=mean_length,sd=mean_length*cv) - #Reset the first length proportion to zero so the first bin includes all - #density smaller than that bin - temp_len_probs[1]<-0 - #subtract the offset length probabilities to calculate the proportion in each - #bin. For each length bin the proportion is how many fish are larger than this - #length but shorter than the next bin length. - temp_len_probs <- c(temp_len_probs[-1],1)-temp_len_probs - length_age_conversion[i,] <- temp_len_probs -} -colnames(length_age_conversion) <- len_bins -rownames(length_age_conversion) <- ages - -#Extract years and fleets from milestone 1 data -start_date <- unique(opaka_dat_fims$timing[opaka_dat_fims$type=="landings"]) -observers <- unique(opaka_dat_fims$name[opaka_dat_fims$type=="length_comp"]) - -#Create data frame for new fleet and year specific length at age conversion proportions -length_age_data <- data.frame( - type = rep("age-to-length-conversion",length(len_bins)*length(ages)*length(observers)*length(start_date)), - name = rep(sort(rep(observers,length(len_bins)*length(ages))),length(start_date)), - age = rep(sort(rep(ages,length(len_bins))),length(observers)*length(start_date)), - length = rep(len_bins,length(ages)*length(observers)*length(start_date)), - timing = rep(start_date,each=length(len_bins)*length(ages)*length(observers)), - value = rep(c(t(length_age_conversion)),length(observers)*length(start_date)), - unit = rep("proportion",length(len_bins)*length(ages)*length(observers)*length(start_date)), - uncertainty = rep(30,length(len_bins)*length(ages)*length(observers)*length(start_date))) - -# Changing the CPUE indices for fleet1 to be a new fleet, fleet4. This helps with model convergence and fitting as it is the longest time series of data available along with the landings. -opaka_dat_fims <- opaka_dat_fims |> - dplyr::mutate(name = ifelse(name == "fleet1" & type == "index", "fleet4", name)) -opaka_dat_fims <- type.convert( - rbind(opaka_dat_fims, length_age_data), - as.is = TRUE -) +opaka_age_to_length <- opaka_dat_fims_raw |> + dplyr::filter(type == "age_to_length_conversion") |> + dplyr::select(-timing) |> + dplyr::mutate( + fleet = "BFISH" + ) |> + tidyr::crossing(timing = alk_years) + +opaka_dat_fims <- opaka_dat_fims_raw |> + dplyr::filter(type != "age_to_length_conversion") |> + dplyr::bind_rows(opaka_age_to_length) |> + dplyr::mutate( + uncertainty = dplyr::if_else( + type == "length_comp" & fleet == "BFISH", + 1, + uncertainty + ) + ) -opaka_dat_fims <- opaka_dat_fims |> - dplyr::mutate(uncertainty = ifelse(type == "length_comp" & name == "fleet2" & value != -999, 1, uncertainty)) |> - dplyr::filter(!(type == "index" & name == "fleet3")) +# TODO clean this up later +# filter to just landings for FRS +# to avoid a bug in FIMS +opaka_dat_fims <- opaka_dat_fims |> + dplyr::filter_out(fleet == "FRS" & type == "index") data_4_model <- FIMS::FIMSFrame(opaka_dat_fims) + ``` The `data_4_model` object contains a `@data` slot that holds a long data frame with: -* 2 fleets: commercial fishery (fleet1) and survey (fleet2) -* landings for fleet 1 -* cpue for fleet 2 -* length composition data for fleet 2 + +* 4 fleets: commercial fishery (FRS) and non-commercial fishery (Non_comm), survey (BFISH), and a new fleet for the CPUE of the commercial fishery (FRS_CPUE) +* landings for FRS and Non_comm +* indices for BFISH and FRS_CPUE +* length composition data for the survey (BFISH) +* weight-at-age data +* age-to-length-conversion data ## Run FIMS model ```{r, max.height='100px', attr.output='.numberLines'} #| label: setup-model - recdevs <- rep$parameters |> - dplyr::filter(stringr::str_detect(Label, "RecrDev")) |> - dplyr::select(Label, Value) + dplyr::filter(stringr::str_detect(Label, "RecrDev")) |> + dplyr::select(Label, Value) -init_naa <- (exp(opaka_ctl$SR_parms["SR_LN(R0)", "INIT"]) * 1000) * exp(-(ages - 1) * 0.135) +init_naa <- (exp(opaka_inputs_bootstrap$ctl$SR_parms["SR_LN(R0)", "INIT"]) * + 1000) * + exp(-(ages - 1) * 0.135) init_naa[n_ages] <- init_naa[n_ages] / 0.135 + # Create default parameters default_parameters <- FIMS::create_default_configurations( - data = data_4_model - ) |> + data = data_4_model +) |> + # # add parametric growth to the configuration (depends on in-development growth module) + # tidyr::unnest(cols = data) |> + # dplyr::mutate( + # module_type = dplyr::if_else( + # module_name == "Growth", + # "VonBertalanffy", + # module_type + # ) + # ) |> + # tidyr::nest(.by = c(model_family, module_name, fleet)) |> FIMS::create_default_parameters( data = data_4_model ) |> - tidyr::unnest(cols = data) |> + tidyr::unnest(cols = data) + +# modify the default parameters +parameters <- default_parameters |> dplyr::rows_update( tibble::tibble( module_name = "Maturity", @@ -365,46 +341,51 @@ default_parameters <- FIMS::create_default_configurations( ), by = c("module_name", "label") ) |> + # # update the growth parameters with initial values from SS3 (depends on in-development growth module) + # dplyr::rows_update( + # tibble::tibble( + # module_name = "Growth", + # label = c( + # "length_at_ref_age_1", + # "length_at_ref_age_2", + # "growth_coefficient_K", + # "reference_age_for_length_1", + # "reference_age_for_length_2" + # ), + # value = c(6, 67.5, .242, 0, 21) + # ), + # by = c("module_name", "label") + # ) |> dplyr::rows_update( tibble::tibble( module_name = "Selectivity", - fleet_name = "fleet1", + fleet = "FRS", label = c("slope", "inflection_point"), # Used age selectivity values value = c(4.5, 1.81), estimation_type = "constant" ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( module_name = "Selectivity", - fleet_name = "fleet2", + fleet = "BFISH", label = c("slope", "inflection_point"), value = c(3, 1), estimation_type = "constant" ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( module_name = "Selectivity", - fleet_name = "fleet3", + fleet = "Non_comm", label = c("slope", "inflection_point"), value = c(4.5, 1.97), estimation_type = "constant" ), - by = c("module_name", "fleet_name", "label") - ) |> - dplyr::rows_update( - tibble::tibble( - module_name = "Selectivity", - fleet_name = "fleet4", - label = c("slope", "inflection_point"), - value = c(4.5, 1.81), - estimation_type = "constant" - ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( @@ -424,17 +405,26 @@ default_parameters <- FIMS::create_default_configurations( ), by = c("module_name", "label", "age") ) |> - # dplyr::rows_update( - # tibble::tibble( - # module_name = "Recruitment", - # # Transformed 0.999 to logit where a previous version just used 0.999 - # # Wondering if we should use logit(0.75) as noted previously for scamp - # # as the null recruitment model? - # label = c("log_rzero", "logit_steep", "log_sd"), - # value = c(opaka_ctl$SR_parms["SR_LN(R0)", "INIT"], -log(1.0 - 0.76) + log(0.76 - 0.2), sca$parm.cons$rec_sigma[8]) - # ), - # by = c("module_name", "label") - # ) |> + dplyr::rows_update( + tibble::tibble( + module_name = "Recruitment", + # Transformed 0.999 to logit where a previous version just used 0.999 + # Wondering if we should use logit(0.75) as noted previously for scamp + # as the null recruitment model? + label = c("log_rzero", "logit_steep", "log_sd"), + value = c( + # log of R0 + log(1000) + # adding log(1000) to account for R0 units in 1000s in SS3 + opaka_inputs_bootstrap$ctl$SR_parms["SR_LN(R0)", "INIT"], + # steepness (with logit transformation) + -log(1.0 - opaka_inputs_bootstrap$ctl$SR_parms["SR_BH_steep", "INIT"]) + + log(opaka_inputs_bootstrap$ctl$SR_parms["SR_BH_steep", "INIT"] - 0.2), + # sigmaR + opaka_inputs_bootstrap$ctl$SR_parms["SR_sigmaR", "INIT"] + ) + ), + by = c("module_name", "label") + ) |> dplyr::rows_update( tibble::tibble( module_name = "Recruitment", @@ -442,51 +432,51 @@ default_parameters <- FIMS::create_default_configurations( time = years[-1], # The last value of the initial numbers at age is the first # recruitment deviation - value = recdevs$Value, + value = recdevs$Value ), by = c("module_name", "label", "time") ) |> dplyr::rows_update( tibble::tibble( module_name = "Fleet", - fleet_name = "fleet1", + fleet = "FRS", time = years, label = "log_Fmort", - value = log(rep$exploitation$FRS) + value = log(rep$exploitation$FRS[-76]) # removing F value for forecast year ), - by = c("module_name", "fleet_name", "label", "time") + by = c("module_name", "fleet", "label", "time") ) |> dplyr::rows_update( tibble::tibble( module_name = "Fleet", - fleet_name = "fleet2", + fleet = "BFISH", label = "log_q", - value = -4.12772 + value = rep$parameters["LnQ_base_BFISH(2)", "Value"] ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( module_name = "Fleet", - fleet_name = "fleet3", + fleet = "Non_comm", label = "log_Fmort", time = years, - value = log(rep$exploitation$Non_comm) - ), - by = c("module_name", "fleet_name", "label", "time") - ) |> - dplyr::rows_update( - tibble::tibble( - module_name = "Fleet", - fleet_name = "fleet4", - label = "log_q", - value = -3.90281 #value from SS + value = log(rep$exploitation$Non_comm[-76]) # removing F value for forecast year ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label", "time") ) -# Run the model with optimization -fit <- default_parameters |> +# # optionally use the defaults instead of the modified parameters +# parameters <- default_parameters +# Run the model without optimization +## Getting error here about mismatch in size of log_sd vector +fit_init <- parameters |> + FIMS::initialize_fims(data = data_4_model) |> + # Model is too big to run on GitHub action if you estimate uncertainty + FIMS::fit_fims(optimize = FALSE, get_sd = FALSE) + +# Run the model with optimization +fit_est <- parameters |> FIMS::initialize_fims(data = data_4_model) |> # Model is too big to run on GitHub action if you estimate uncertainty FIMS::fit_fims(optimize = TRUE, get_sd = FALSE) @@ -494,117 +484,7 @@ fit <- default_parameters |> ## Plotting Results -```{r} -#| label: comparison-plots - -index_results <- data.frame( - observed = FIMS::m_index(data_4_model, "fleet2"), - expected = FIMS::get_report(fit)[["index_expected"]][[2]] -) |> -dplyr::mutate(year = years) |> -dplyr::filter(year > 2016) -#print(index_results) - -ggplot2::ggplot(index_results, ggplot2::aes(x = year, y = observed)) + - ggplot2::geom_point() + - ggplot2::xlab("Year") + - ggplot2::ylab("Index (mt)") + - ggplot2::geom_line(ggplot2::aes(x = year, y = expected), color = "blue") + - ggplot2::theme_bw() - -cpue_results <- data.frame( - observed = FIMS::m_index(data_4_model, "fleet4"), - expected = FIMS::get_report(fit)[["index_expected"]][[2]] -) |> -dplyr::mutate(year = years) |> -dplyr::filter(observed >0) -#print(cpue_results) -ggplot2::ggplot(cpue_results, ggplot2::aes(x = year, y = observed)) + - ggplot2::geom_point() + - ggplot2::xlab("Year") + - ggplot2::ylab("CPUE") + - ggplot2::geom_line(ggplot2::aes(x = year, y = expected), color = "blue") + - ggplot2::theme_bw() - - -catch_results <- data.frame( - observed = c(FIMS::m_landings(data_4_model, fleet = "fleet1"), FIMS::m_landings(data_4_model, fleet = "fleet3")), - expected = c(FIMS::get_report(fit)[["landings_expected"]][[1]], FIMS::get_report(fit)[["landings_expected"]][[3]]), - fleet = rep(c("fleet1", "fleet3"), each = 75) -) |> -dplyr::mutate(year = rep(years, 2)) -#print(catch_results) - -ggplot2::ggplot(catch_results, ggplot2::aes(x = year, y = observed)) + - ggplot2::geom_point(ggplot2::aes(color = fleet)) + - ggplot2::xlab("Year") + - ggplot2::ylab("Catch (mt)") + - ggplot2::geom_line(ggplot2::aes(x = year, y = expected, color = fleet)) + - ggplot2::theme_bw() - -biomass <- rep$timeseries |> -dplyr::select(Yr, SpawnBio, Bio_all) |> -dplyr::filter(Yr > 1947) |> ##CHECK: including "initial year" to match length with FIMS but need to check on FIMS -dplyr::rename("SS_SpawnBio" = "SpawnBio", - "SS_Bio" = "Bio_all") |> -dplyr::mutate(FIMS_SpawnBio = FIMS::get_report(fit)[["spawning_biomass"]][[1]] , - FIMS_Bio = FIMS::get_report(fit)[["biomass"]][[1]]) |> ##CHECK: Is FIMS ssb reporting n_years+1 or initial year-1? -tidyr::pivot_longer(cols = -Yr) |> -tidyr::separate_wider_delim(cols = "name", delim = "_", names = c("Model", "Type")) - -ggplot2::ggplot(biomass, ggplot2::aes(x = Yr, y = value)) + - ggplot2::geom_line(ggplot2::aes(color = Model)) + - ggplot2::xlab("Year") + - ggplot2::ylab("") + - ggplot2::facet_wrap(~Type, scales = "free_y") + - ggplot2::theme_bw() - -recruits <- rep$recruit |> - dplyr::select(Yr, exp_recr, raw_dev) |> - dplyr::rename("SS_recruit" = "exp_recr", - "SS_recdev" = "raw_dev", - "Year" = "Yr") |> - dplyr::mutate(FIMS_recruit = FIMS::get_report(fit)[["expected_recruitment"]][[1]][1:75], - FIMS_recdev = c( - NA, - FIMS::get_estimates(fit) |> - dplyr::filter(label == "log_devs", module_name == "Recruitment") |> - dplyr::pull(estimated) - ), - SS_recruit = SS_recruit * 1000) |> - tidyr::pivot_longer(cols = -Year) |> - tidyr::separate_wider_delim(cols = "name", delim = "_", names = c("Model", "Type")) - -ggplot2::ggplot(recruits, ggplot2::aes(x = Year, y = value, color = Model)) + - ggplot2::geom_line() + - ggplot2::facet_wrap(~Type, scales = "free_y") + - ggplot2::theme_bw() -``` - - -```{r} -#| eval: false -#| label: proportions-plots-not-run -## Checking fit to proportion catch number at length -pcnal <- matrix(data = fit@report$pcnal[[2]], nrow = nlengths) - -prop.dat <- opaka_dat_fims |> -filter(type == "length" & name == "fleet2") |> -group_by(datestart) |> -reframe(prop = value/sum(value)) -pcnal.obs <- matrix(data = prop.dat$prop, nrow = nlengths) -head(pcnal.obs) - -plot(x = 1:nlengths, y = pcnal.obs[,73], pch = 16, ylim = c(0,1)) -lines(x = 1:nlengths, y = pcnal[,73]) -pcnal[,69] - -## checking estimated numbers at age -head(fit@report$naa[[1]]) -naa_mat <- matrix(data = fit@report$naa[[1]], nrow = n_ages) -head(naa_mat) - -``` +TODO: add plots in the future ```{r} #| label: cleanup diff --git a/content/R/NWFSC-petrale_simplify_assessment.R b/content/R/NWFSC-petrale_simplify_assessment.R new file mode 100644 index 0000000..0b72988 --- /dev/null +++ b/content/R/NWFSC-petrale_simplify_assessment.R @@ -0,0 +1,282 @@ +#' Simplify the NWFSC petrale sole SS3 assessment inputs +#' +#' Simplifies the original petrale sole production assessment model to better +#' match the options currently available in FIMS. This function is not intended +#' to be run by FIMS users, rather to document the process of simplifying the +#' SS3 model for future reference. As the set of FIMS features expands, this +#' script could be revised to do fewer simplification steps. +#' +#' @param new_model_dir Directory to write the simplified SS3 input files to. +#' @return Invisibly returns the modified SS3 input list after writing the +#' simplified files to `new_model_dir`. + +NWFSC_petrale_simplify_assessment <- function( + new_model_dir = "models/2023.a050.003_FIMS_case-study_wtatage" +) { + # read SS3 input files from petrale sole assessment on GitHub + petrale_input <- r4ss::SS_read( + "https://raw.githubusercontent.com/pfmc-assessments/petrale/main/models/2023.a034.001/", + read_wtatage = TRUE + ) + + ################################################################## + ## SIMPLIFY SS3 DATA FILE + + dat <- petrale_input$dat + + # change dimensions + dat$Nsexes <- 1 # previously 2 + dat$Nages <- 17 # previously 40 + + # filter indices (just include WCGBTS) + dat$CPUE <- dat$CPUE |> + dplyr::filter(index == 4) + + # remove discard info + dat$N_discard_fleets <- 0 + dat$discard_fleet_info <- NULL + dat$discard_data <- NULL + + # remove mean body weight + dat$use_meanbodywt <- 0 + dat$meanbodywt <- NULL + dat$DF_for_meanbodywt <- NULL + + # simplify ageing error matrices to just use first matrix + # and only for new range of ages + dat$N_ageerror_definitions <- 1 + dat$ageerror <- dat$ageerror[1:2, paste0("age", 0:dat$Nages)] # first row is mean, second row is sd + + # simplify age comps + dat$agecomp <- + dat$agecomp |> + # exclude fleet 3, and only use marginal data for fleet 4 + dplyr::filter(fleet %in% c(1:2, -4)) |> + # exclude male compositions (columns matching pattern "m" followed + # by numeric value for length bin, but doesn't match "month") + dplyr::select(!dplyr::matches("^m\\d+$")) |> + # assign observations to sex 1, ageing error 1, and include previously excluded marginal ages + dplyr::mutate(sex = 1, ageerr = 1, fleet = abs(fleet)) |> + # remove redundant observations + # NOTE: this removes some age comp data because there + # were years with multiple observations from the same fleet + # due to multiple ageing error matrices + dplyr::distinct(year, fleet, .keep_all = TRUE) + + # simplify length comps + dat$lencomp <- + dat$lencomp |> + # exclude fleet 3 and any discard length comps (part == 1) + dplyr::filter(part %in% c(0, 2) & fleet %in% c(1, 2, 4)) |> + # exclude male compositions (columns matching pattern "m" followed + # by numeric value for length bin, but doesn't match "month") + dplyr::select(!dplyr::matches("^m\\d+$")) |> + dplyr::mutate(sex = 1) + + ################################################################## + ## SIMPLIFY SS3 CONTROL FILE + + ctl <- petrale_input$ctl + ctl$Nages <- dat$Nages # previously 40 + ctl$Nsexes <- dat$Nsexes # previously 2 + + # remove all male mortality and growth parameters + ctl$MG_parms <- ctl$MG_parms |> + dplyr::filter(!grepl("_Mal_", rownames(ctl$MG_parms))) + ctl$max_bias_adj <- -1 # set bias adjust = 1.0 for all years + + # fix growth parameters by setting phase negative for everything except M + ctl$MG_parms <- ctl$MG_parms |> + dplyr::mutate( + PHASE = ifelse(grepl("NatM_", rownames(ctl$MG_parms)), PHASE, -abs(PHASE)) + ) + + # # fix M at 0.1 in case this helps with estimation problems + # ctl$MG_parms <- ctl$MG_parms |> + # dplyr::mutate( + # PHASE = ifelse( + # grepl("NatM_", rownames(ctl$MG_parms)), + # -abs(PHASE), + # PHASE + # ), + # INIT = ifelse(grepl("NatM_", rownames(ctl$MG_parms)), 0.1, INIT) + # ) + + # filter inputs stuff to only fleet 4 (WCGBTS) + ctl$Q_options <- ctl$Q_options |> + dplyr::filter(fleet == 4) + ctl$Q_parms <- ctl$Q_parms |> + dplyr::filter(grepl("WCGBTS", rownames(ctl$Q_parms))) + + # add new age-based selectivity parameters + # age-based logistic (pattern 12) for fleets 1 and 4 + # have south fishery (fleet 2) mirror fleet 1 + # no parameters for fleet 3 (ignored index left in place to avoid renumbering fleets) + ctl$age_selex_types <- ctl$age_selex_types |> + dplyr::mutate( + Pattern = c(12, 15, 0, 12), + Special = c(0, 1, 0, 0) + ) # fleet 2 mirror fleet 1 + + # create 4 parameter rows + ctl$age_selex_parms <- cbind( + data.frame( + LO = 1, + HI = 10, + INIT = c(5, 2, 5, 2), + PRIOR = 0, + PR_SD = 99, + PR_type = 0, + PHASE = 2 + ), + matrix(0, nrow = 4, ncol = 7) # this fills in columns 8 to 14 with zeros + ) + names(ctl$age_selex_parms)[8:14] <- names(ctl$size_selex_parms)[8:14] + rownames(ctl$age_selex_parms) <- paste0( + "AgeSel_", + c("P_1_", "P_2_"), + c("fleet1", "fleet1", "fleet4", "fleet4") + ) + + # remove length-based selectivity + ctl$size_selex_types <- ctl$size_selex_types |> + dplyr::mutate(Pattern = 0, Discard = 0, Male = 0) + ctl$size_selex_parms <- NULL + ctl$size_selex_parms_tv <- NULL + + # remove blocks + ctl$N_Block_Designs <- 0 + ctl$blocks_per_pattern <- NULL + ctl$Block_Design <- NULL + + # turn off variance adjustments + ctl$DoVar_adjust <- 0 + ctl$Variance_adjustment_list <- NULL + # use wtatage file instead of growth parameters + ctl$EmpiricalWAA <- 1 + + # put all recdevs into the main vector + ctl$MainRdevYrFirst <- 1861 + ctl$MainRdevYrLast <- 2022 + ctl$recdev_early_start <- 1845 + ctl$recdev_early_phase <- -3 + ctl$Fcast_recr_phase <- -4 + + ctl$recdev_phase <- 5 + + ################################################################## + ## SIMPLIFY wtatage file + + wtatage <- petrale_input$wtatage + + if (FALSE) { + # explore estimated wtatage from SS3 model (differs among fleets due to + # length-based selectivity and changes over time due to selectivity blocks) + # nevertheless, it's easiest to just use the population wtatage for all fleets + library(ggplot2) + wtatage |> + dplyr::filter(fleet %in% -1:4, sex == 1) |> + select(year, fleet, "20") |> + rename(wt = "20") |> + ggplot(aes(x = year, y = wt, color = as.factor(fleet))) + + geom_line() + } + + # wtatage for fleets -1:0 (population), and fleets 1:4 (real fleets) + wtatage_pop <- wtatage |> + dplyr::select(1:6, paste(0:dat$Nages)) |> + dplyr::filter(year == 1876 & sex == 1 & fleet == -1) |> + dplyr::mutate(year = -year) + # original model had a fecundity relationship + # so recalculate the maturity * fecundity based on the population weight-at-age + # and an age-based maturity curve (logistic with 50% maturity at age 6 and slope of 1.5) + # this is a rough approximation to the length-based maturity converted to age + # within the original model, but should be sufficient for testing the simplified model + maturity <- 1 / (1 + exp(-1.5 * (0:dat$Nages - 6))) + matfec <- maturity * + as.numeric(dplyr::select(wtatage_pop, paste0(0:dat$Nages))) + wtatage_matfec <- wtatage_pop + wtatage_matfec[, paste(0:dat$Nages)] <- matfec + + # # old code to use SS3 output for maturity * fecundity + # wtatage_matfec <- wtatage |> + # dplyr::select(1:6, paste(0:dat$Nages)) |> + # dplyr::filter(year == 1876 & sex == 1 & fleet == -2) |> + # dplyr::mutate(year = -year) + wtatage_simple <- rbind( + wtatage_matfec, # fleet = -2 + wtatage_pop, # fleet = -1 + wtatage_pop, # fleet = 0 + wtatage_pop, # fleet = 1 + wtatage_pop, # fleet = 2 + wtatage_pop, # fleet = 3 + wtatage_pop # fleet = 4 + ) |> + dplyr::mutate(fleet = -2:4) + + # simplify forecast file + fore <- petrale_input$fore + fore$Forecast <- 0 # turn off forecast + fore$Bmark_relF_Basis <- 1 # use year range for benchmark relative F + + # write modified data and control files + new_input <- petrale_input + new_input$dat <- dat + new_input$ctl <- ctl + new_input$wtatage <- wtatage_simple + new_input$fore <- fore + + # reduce run display detail + new_input$start$run_display_detail <- 0 + + new_input$dir <- new_model_dir + # create directory if it doesn't exist + if (!dir.exists(new_input$dir)) { + dir.create(new_input$dir) + } + + r4ss::SS_write(new_input, dir = new_input$dir, overwrite = TRUE) + + # run model and compare output + if (FALSE) { + r4ss::run( + new_input$dir, + show_in_console = TRUE, + skipfinished = FALSE, + extras = "-maxfn 3000" # was exceeding default maxfn 100 in early phases + ) + r4ss::run( + new_input$dir, + show_in_console = TRUE, + skipfinished = FALSE, + extras = "-hess_step" # was exceeding default maxfn 100 in early phases + ) + + # read output from model that was just run + petrale_simple_output <- r4ss::SS_output( + new_input$dir, + printstats = FALSE, + verbose = FALSE, + covar = FALSE + ) + + r4ss::SS_plots(petrale_simple_output) + + # read SS3 output files from petrale sole assessment on GitHub + petrale_output1 <- r4ss::SS_output( + "https://raw.githubusercontent.com/pfmc-assessments/petrale/main/models/2023.a034.001", + printstats = FALSE, + verbose = FALSE, + covar = FALSE + ) + + # compare outputs from simplified and original SS3 models + # (fairly large changes in spawning output, perhaps due to removal of the older age classes) + SSplotComparisons(SSsummarize(list( + petrale_simple_output, + petrale_output1 + ))) + } + + invisible(new_input) +} diff --git a/content/R/get_asap_data.R b/content/R/get_asap_data.R index d865279..7c3ddfa 100644 --- a/content/R/get_asap_data.R +++ b/content/R/get_asap_data.R @@ -1,7 +1,7 @@ # need to think about how to deal with multiple fleets, only using 1 right now get_asap_data <- function(asap_input) { res <- data.frame(type = character(), - name = character(), + fleet = character(), age = integer(), timing = double(), value = double(), @@ -10,7 +10,7 @@ get_asap_data <- function(asap_input) { years_in_model <- seq(asap_input$parms$styr, asap_input$parms$endyr) landings <- data.frame(type = "landings", - name = "fishery", + fleet = "fishery", age = NA, timing = years_in_model, value = as.numeric(asap_input$catch.obs[1,]), @@ -26,7 +26,7 @@ get_asap_data <- function(asap_input) { for (i in seq(asap_input$parms$nindices)) { index <- data.frame( type = "index", - name = paste0("survey", i), + fleet = paste0("survey", i), age = NA_integer_, timing = years_in_model[asap_input$index.year.counter[[i]]], value = as.numeric(asap_input$index.obs[[i]]), @@ -46,7 +46,7 @@ get_asap_data <- function(asap_input) { catchage <- data.frame( type = "age_comp", - name = "fishery", + fleet = "fishery", age = rep(seq(1,asap_input$parms$nages), asap_input$parms$nyears), timing = rep(seq(asap_input$parms$styr, asap_input$parms$endyr), each = asap_input$parms$nages), value = as.numeric(t(asap_input$catch.comp.mats$catch.fleet1.ob)), @@ -57,7 +57,7 @@ get_asap_data <- function(asap_input) { # loop over all indices for (i in 1:asap_input$parms$nindices){ indexage <- data.frame(type = "age_comp", - name = paste0("survey", i), + fleet = paste0("survey", i), age = rep(seq(1,asap_input$parms$nages), asap_input$parms$nyears), timing = rep(seq(asap_input$parms$styr, asap_input$parms$endyr), each = asap_input$parms$nages), value = as.numeric(t(asap_input$index.comp.mats[[i*2-1]])), @@ -84,7 +84,7 @@ get_asap_data <- function(asap_input) { weight_at_age <- data.frame( type = "weight_at_age", - name = NA_character_, + fleet = NA_character_, age = rep( seq(asap_input[["parms"]][["nages"]]), length(seq(asap_input[["parms"]][["styr"]], asap_input[["parms"]][["endyr"]] + 1)) diff --git a/content/R/get_ss3_data.R b/content/R/get_ss3_data.R deleted file mode 100644 index 13f2c53..0000000 --- a/content/R/get_ss3_data.R +++ /dev/null @@ -1,220 +0,0 @@ -#' Convert SS3 data into format required by FIMS (works for petrale and opaka so far) -#' -#' Uses output from `r4ss::SS_read()` or `r4ss::SS_readdat()` and does -#' filtering, simplifying, and reformatting. -#' -#' @param ss3_inputs A list containing `dat` and `wtatage` such as that -#' created by `r4ss::SS_read()` or by running `r4ss::SS_readdat()` and -#' `r4ss::SS_readwtatage()` and combining the results in a list. -#' Note: if the SS3 model has parametric growth then `r4ss::SS_read()` won't -#' read the `wtatage.ss` file and it needs to be added to the list by -#' separately running `r4ss::SS_readwtatage()` or taking it from the list -#' created by `r4ss::SS_output()`. -#' -#' @param fleets Which fleets to include in the processed output. -#' Note that the only start year population weight-at-age is read from the -#' `wtatage` element (fleet = 0). NULL will default to including all fleets -#' from the SS3 model. -#' @param ages Vector of ages to index. NULL will default to using -#' all age data bins from the SS3 model. -#' @param lengths Vector of lengths to index. NULL will default to using -#' all length data bins from the SS3 model. -#' @return A data frame that can be passed to `FIMS::FIMSFrame()` -#' @author Ian G. Taylor, Megumi Oshima, Kelli F. Johnson -#' @export - -get_ss3_data <- function(ss3_inputs, fleets = NULL, ages = NULL, lengths = NULL) { - # check inputs for necessary elements - if (!is.list(ss3_inputs) || !"dat" %in% names(ss3_inputs)) { - stop("`ss3_inputs` should be a list containing both 'dat' and 'wtatage'") - } - if (!"wtatage" %in% names(ss3_inputs)) { - stop("'ss3_inputs' is missing element 'wtatage'. You may have to add it by running 'r4ss::SS_readwtatage()'") - } - if (any(ss3_inputs[["wtatage"]][["year"]] < 0)) { - stop( - "The 'wtatage' element includes negative years which can't yet be processed by this function, please use the wtatage.ss_new file." - ) - } - - # pull out dat element from the list to simplify code - dat <- ss3_inputs$dat - - # fill in any missing inputs - if (is.null(fleets)) { - fleets <- seq_along(dat$fleetnames) - } - if (is.null(ages)) { - ages <- dat$agebin_vector - } - if (is.null(lengths)) { - lengths <- dat$lbin_vector - } - - # create empty data frame - res <- data.frame( - type = character(), - name = character(), - age = integer(), - length = integer(), - timing = character(), - value = double(), - unit = character(), - uncertainty = double() - ) - - # aggregate landings across fleets - catch_by_year_fleet <- dat$catch |> - dplyr::filter(year != -999) |> # year = -999 in SS3 designates initial equilibrium catch - dplyr::filter(fleet %in% fleets) - - # convert landings to FIMSFrame format - landings <- data.frame( - type = "landings", - name = paste0("fleet", catch_by_year_fleet$fleet), # landings aggregated to fleet 1 - age = NA, - length = NA, - timing = catch_by_year_fleet$year, - value = catch_by_year_fleet$catch, - unit = "mt", - uncertainty = catch_by_year_fleet$catch_se - ) - - # check for any gaps in landings time series - years <- min(catch_by_year_fleet$year):max(catch_by_year_fleet$year) - if (!all(years %in% catch_by_year_fleet$year)) { - stop("missing years in landings") - } - - # convert indices to FIMSFrame format - index_info <- dat$CPUE |> - dplyr::filter(index %in% fleets) |> - dplyr::select(year, index, obs, se_log) |> - dplyr::arrange(index, year) - - indices <- data.frame( - type = "index", - name = paste0("fleet", index_info$index), - age = NA, - length = NA, - timing = index_info$year, - value = index_info$obs, - unit = "", - uncertainty = index_info$se_log - ) - - if (!is.null(dat$agecomp)) { - # partially convert age comps (filter, make into long table) - - # first rescale females to sum to 1.0 - # (data processing step had females + males sum to 100 for no good reason) - dat$agecomp$sum_fem <- - dat$agecomp |> - dplyr::select(dplyr::starts_with(c("f", "a"), ignore.case = FALSE)) |> # get female comps (or comps if single-sex) - rowSums() - # couldn't figure out dplyr approach to rescaling the subset of columns - # with female proportions to sum to 1.0 - fcols <- dat$agecomp |> dplyr::select(dplyr::starts_with("f", ignore.case = FALSE)) - if (length(fcols) > 0) { - dat$agecomp[, names(dat$agecomp) %in% paste0("f", ages)] <- - dat$agecomp[, names(dat$agecomp) %in% paste0("f", ages)] / - dat$agecomp$sum_fem - } else { - dat$agecomp[, names(dat$agecomp) %in% paste0("a", ages)] <- - dat$agecomp[, names(dat$agecomp) %in% paste0("a", ages)] / - dat$agecomp$sum_fem - } - - # further processing - age_info <- - dat$agecomp |> - dplyr::filter(fleet %in% fleets) |> # filter by requested fleets - dplyr::mutate(fleet = abs(fleet)) |> # convert any negative fleet to positive - dplyr::select(!dplyr::matches("^m[0-9]")) |> # exclude male comps - tidyr::pivot_longer( # convert columns f1...f17 to values in a new "age" colum of a longer table - cols = dplyr::matches("^f[0-9]") | dplyr::matches("^a[0-9]"), # 2-sex model uses f1, f2, ...; 1-sex model uses a1, a2, ... - names_to = "age", - values_to = "value" - ) |> - dplyr::mutate(age = as.numeric(substring(age, first = 2))) |> # convert "f17" to 17 - dplyr::select(year, fleet, Nsamp, age, value) |> - # Find missing age in composition data and fill in for each fleet - dplyr::group_by(fleet, year) |> - tidyr::complete(age = ages) |> - tidyr::fill(Nsamp, .direction = "updown") |> - dplyr::mutate(value = ifelse(is.na(value), 0, value)) |> - dplyr::ungroup() |> - dplyr::arrange(fleet, year, age) - - # finish converting age comps to FIMSFrame format - agecomps <- data.frame( - type = "age_comp", - name = paste0("fleet", abs(age_info$fleet)), # abs to include fleet == -4 - age = age_info$age, - length = NA, - timing = age_info$year, - value = age_info$value + 0.001, # add constant to avoid 0 values - unit = "", - # Q: should uncertainty here be the total sample size across bins, or the samples within the bin? - # uncertainty = round(age_info$Nsamp * age_info$value) - uncertainty = round(age_info$Nsamp) - ) - } else { - agecomps <- NULL - } - - ## Length composition data - if (!is.null(dat[["lencomp"]])) { - # leaving out the re-scaling part for females to 1 - len_info <- - dat$lencomp |> - dplyr::filter(fleet %in% fleets) |> # filter by requested fleets - dplyr::mutate(fleet = abs(fleet)) |> # convert any negative fleet to positive - dplyr::select(!dplyr::matches("^m[0-9]")) |> # exclude male comps - tidyr::pivot_longer( # convert columns f1...f17 to values in a new "length" colum of a longer table - cols = dplyr::matches("^f[0-9]") | dplyr::matches("^l[0-9]"), # 2-sex model uses f1, f2, ...; 1-sex model uses a1, a2, ... - names_to = "length", - values_to = "value" - ) |> - dplyr::mutate(length = as.numeric(substring(length, first = 2))) |> # convert "l17" to 17 - dplyr::select(year, fleet, Nsamp, length, value) |> - dplyr::arrange(fleet, year, length) - - # finish converting age comps to FIMSFrame format - lencomps <- data.frame( - type = "length_comp", # will likely need to change name - name = paste0("fleet", abs(len_info$fleet)), # abs to include fleet == -4 - age = NA, - length = len_info$length, - timing = len_info$year, - value = len_info$value + 0.001, # add constant to avoid 0 values - unit = "", - # Q: should uncertainty here be the total sample size across bins, or the samples within the bin? - # uncertainty = round(len_info$Nsamp * len_info$value) - uncertainty = round(len_info$Nsamp) - ) - } else { - lencomps <- NULL # not sure if we need this but wanting to avoid an error if missing age or length comps - } - - ## Weight-at-age data - wtatage <- ss3_inputs$wtatage |> - dplyr::filter(fleet == 0 & sex == 1 & seas == 1 & birthseas == 1) |> - dplyr::select("year", dplyr::matches("[0-9]+")) |> - tidyr::pivot_longer(names_to = "age", cols = -year) |> - dplyr::filter(age %in% ages) |> - dplyr::mutate( - type = "weight-at-age", - name = "fleet1", - age = as.integer(age), - length = NA, - timing = year, - value = value / 1000, # covert to metric tons (SS3) - unit = "mt", - uncertainty = NA - ) |> - dplyr::select(-year) - - # combine all data sources - res <- rbind(res, landings, indices, agecomps, lencomps, wtatage) -} diff --git a/content/R/get_ss3_timeseries.R b/content/R/get_ss3_timeseries.R deleted file mode 100644 index 1c2ef45..0000000 --- a/content/R/get_ss3_timeseries.R +++ /dev/null @@ -1,48 +0,0 @@ -#' Extract time series output from SS3 model and put into a long data frame -#' -#' Simple function to reformat output from the timeseries table returned by -#' `r4ss::SS_output()` to facilitate use of tidyverse style functions and -#' facilitate comparison with FIMS output. Currently only extracts F for a -#' single fleet. -#' -#' @param model The output from `r4ss::SS_output()` -#' @param platform A user-specified label to differentiate from additional -#' output associated with FIMS or other platforms -#' @return A long data frame (actually a tibble) with time series output for -#' 4 quantities (so far) -#' @author Ian G. Taylor -#' @examples -#' \dontrun{ -#' # read SS3 models from location on Ian's computer -#' p1 <- r4ss::SS_output("c:/ss/Petrale/Petrale2023/petrale/models/2023.a034.001/") -#' p2 <- r4ss::SS_output("c:/ss/Petrale/Petrale2023/petrale/models/2023.a050.002_FIMS_case-study_wtatage/") -#' # # saving all model output creates large files (7MB for original) -#' # saveRDS(p1, file = "content/data_files/NWFSC-petrale-SS3-original.rds") -#' # saveRDS(p2, file = "content/data_files/NWFSC-petrale-SS3-simplified.rds") -#' # combine SS3 model time series into data frame using function above -#' timeseries_compare <- rbind( -#' get_ss3_timeseries(model = p1, platform = "ss3_original"), -#' get_ss3_timeseries(model = p2, platform = "ss3_simplified") -#' ) -#' # save data frame of time series results to compare with FIMS -#' saveRDS(timeseries_compare, file = "content/data_files/NWFSC-petrale-SS3-timeseries.rds") -#' } - -get_ss3_timeseries <- function(model, platform = "ss3") { - timeseries_ss3 <- model$timeseries |> - dplyr::filter(Yr %in% timeseries$year) |> # filter for matching years only (no forecast) - dplyr::select(Yr, Bio_all, SpawnBio, Recruit_0, "F:_1") |> # select quants of interest - dplyr::rename( # change to names used with FIMS - year = Yr, - biomass = Bio_all, ssb = SpawnBio, recruitment = Recruit_0, F_mort = "F:_1" - ) |> - dplyr::mutate(ssb = 1000 * ssb) |> - tidyr::pivot_longer( # convert quantities in separate columns into a single value column - cols = -1, - names_to = "type", - values_to = "value" - ) |> - dplyr::arrange(type) |> # sort by type instead of year - dplyr::mutate(platform = platform) - return(timeseries_ss3) -} diff --git a/content/R/pk_prepare_FIMS_inputs.R b/content/R/pk_prepare_FIMS_inputs.R index c40923c..b306d45 100644 --- a/content/R/pk_prepare_FIMS_inputs.R +++ b/content/R/pk_prepare_FIMS_inputs.R @@ -164,7 +164,7 @@ prepare_pollock_data <- function( ## put into fims friendly form res <- data.frame( type = character(), - name = character(), + fleet = character(), age = integer(), timing = double(), value = double(), @@ -173,7 +173,7 @@ prepare_pollock_data <- function( ) landings <- data.frame( type = "landings", - name = "fleet1", + fleet = "fleet1", age = NA, timing = seq(fimsdat$styr, fimsdat$endyr), value = as.numeric(fimsdat$cattot) * 1e3, @@ -182,7 +182,7 @@ prepare_pollock_data <- function( ) index2 <- data.frame( type = "index", - name = "survey2", + fleet = "survey2", age = NA, timing = seq(fimsdat$styr, fimsdat$endyr), value = ifelse(ind2 > 0, ind2 * 1e9, ind2), @@ -191,7 +191,7 @@ prepare_pollock_data <- function( ) index3 <- data.frame( type = "index", - name = "survey3", + fleet = "survey3", age = NA, timing = seq(fimsdat$styr, fimsdat$endyr), value = ifelse(ind3 > 0, ind3 * 1e9, ind3), @@ -200,7 +200,7 @@ prepare_pollock_data <- function( ) index6 <- data.frame( type = "index", - name = "survey6", + fleet = "survey6", age = NA, timing = seq(fimsdat$styr, fimsdat$endyr), value = ifelse(ind6 > 0, ind6 * 1e9, ind6), @@ -210,7 +210,7 @@ prepare_pollock_data <- function( ## these have -999 for missing data years catchage <- data.frame( type = "age_comp", - name = "fleet1", + fleet = "fleet1", age = rep(seq(1, n_ages), n_years), timing = rep( seq(fimsdat$styr, fimsdat$endyr), @@ -222,7 +222,7 @@ prepare_pollock_data <- function( ) indexage2 <- data.frame( type = "age_comp", - name = "survey2", + fleet = "survey2", age = rep(seq(1, n_ages), n_years), timing = rep( seq(fimsdat$styr, fimsdat$endyr), @@ -234,7 +234,7 @@ prepare_pollock_data <- function( ) indexage3 <- data.frame( type = "age_comp", - name = "survey3", + fleet = "survey3", age = rep(seq(1, n_ages), n_years), timing = rep( seq(fimsdat$styr, fimsdat$endyr), @@ -246,7 +246,7 @@ prepare_pollock_data <- function( ) indexage6 <- data.frame( type = "age_comp", - name = "survey6", + fleet = "survey6", age = rep(seq(1, n_ages), n_years), timing = rep( seq(fimsdat$styr, fimsdat$endyr), @@ -269,10 +269,11 @@ prepare_pollock_data <- function( each = n_ages ) ) - weightsfishery <- rbind( - pkinput$dat$wt_srv1, - pkinput$dat$wt_srv1[NROW(pkinput$dat$wt_srv1), ] - ) + weightsfishery <- do.call(rbind, replicate( + length(years) + 1, + pkinput$dat$wt_srv1[1, ], + simplify = FALSE + )) colnames(weightsfishery) <- ages rownames(weightsfishery) <- c(years, max(years) + 1) weightatage_data <- tidyr::pivot_longer( diff --git a/content/SEFSC-scamp.qmd b/content/SEFSC-scamp.qmd index fc1bdd1..8d8586c 100644 --- a/content/SEFSC-scamp.qmd +++ b/content/SEFSC-scamp.qmd @@ -114,7 +114,7 @@ survey_ac[!is.na(sca$t.series$acomp.CVT.n), ] <- sca$comp.mats$acomp.CVT.ob ## put data into fims friendly form fleet1_landings_df <- data.frame( type = "landings", - name = "fleet1", + fleet = "fleet1", age = NA, timing = seq(styr, endyr), value = as.numeric(fleet1_landings), @@ -124,7 +124,7 @@ fleet1_landings_df <- data.frame( fleet2_landings_df <- data.frame( type = "landings", - name = "fleet2", + fleet = "fleet2", age = NA, timing = seq(styr, endyr), value = as.numeric(fleet2_landings), @@ -134,7 +134,7 @@ fleet2_landings_df <- data.frame( survey_index_df <- data.frame( type = "index", - name = "survey1", + fleet = "survey1", age = NA, timing = seq(styr, endyr), value = as.numeric(survey_index), @@ -144,7 +144,7 @@ survey_index_df <- data.frame( fleet1_ac_df <- data.frame( type = "age_comp", - name = "fleet1", + fleet = "fleet1", age = rep(seq(1, n_ages), n_years), timing = rep(seq(styr, endyr), each = n_ages), value = as.numeric(t(fleet1_ac)), @@ -154,7 +154,7 @@ fleet1_ac_df <- data.frame( fleet2_ac_df <- data.frame( type = "age_comp", - name = "fleet2", + fleet = "fleet2", age = rep(seq(1, n_ages), n_years), timing = rep(seq(styr, endyr), each = n_ages), value = as.numeric(t(fleet2_ac)), @@ -164,7 +164,7 @@ fleet2_ac_df <- data.frame( survey_ac_df <- data.frame( type = "age_comp", - name = "survey1", + fleet = "survey1", age = rep(seq(1, n_ages), n_years), timing = rep(seq(styr, endyr), each = n_ages), value = as.numeric(t(survey_ac)), @@ -174,7 +174,7 @@ survey_ac_df <- data.frame( weight_at_age <- data.frame( type = "weight_at_age", - name = NA_character_, + fleet = NA_character_, age = seq(n_ages), timing = NA_integer_, value = sca$a.series$wgt.mt, @@ -220,32 +220,32 @@ default_parameters <- FIMS::create_default_configurations(data_4_model) |> dplyr::rows_update( tibble::tibble( module_name = "Selectivity", - fleet_name = "fleet1", + fleet = "fleet1", label = c("slope", "inflection_point"), value = c(sca$parm.cons$selpar_slope_COM2[8], sca$parm.cons$selpar_A50_COM2[8]), estimation_type = "fixed_effects" ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( module_name = "Selectivity", - fleet_name = "fleet2", + fleet = "fleet2", label = c("slope", "inflection_point"), value = c(sca$parm.cons$selpar_slope1_REC2[8], sca$parm.cons$selpar_A50_REC2[8]), estimation_type = "fixed_effects" ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( module_name = "Selectivity", - fleet_name = "survey1", + fleet = "survey1", label = c("slope", "inflection_point"), value = c(sca$parm.cons$selpar_slope1_CVT[8], sca$parm.cons$selpar_A501_CVT[8]), estimation_type = "fixed_effects" ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( @@ -295,34 +295,34 @@ default_parameters <- FIMS::create_default_configurations(data_4_model) |> # dplyr::rows_update( # tibble::tibble( # module_name = "Fleet", - # fleet_name = "fleet1", + # fleet = "fleet1", # time = years, # label = "log_Fmort", # value = log(sca$t.series$F.COM) # ), - # by = c("module_name", "fleet_name", "label", "time") + # by = c("module_name", "fleet", "label", "time") # ) |> # dplyr::rows_update( # tibble::tibble( # module_name = "Fleet", - # fleet_name = "fleet2", + # fleet = "fleet2", # time = years, # label = "log_Fmort", # value = log(sca$t.series$F.REC) # ), - # by = c("module_name", "fleet_name", "label", "time") + # by = c("module_name", "fleet", "label", "time") # ) |> dplyr::rows_update( tibble::tibble( module_name = "Fleet", - fleet_name = "survey1", + fleet = "survey1", # time = years, time = seq(get_start_year(data_4_model), get_end_year(data_4_model)), label = "log_Fmort", # value = log(sca$parms$q.CVT) value = -200 ), - by = c("module_name", "fleet_name", "label", "time") + by = c("module_name", "fleet", "label", "time") ) #################################################################################### @@ -383,7 +383,7 @@ index_results <- data.frame( observed = c( rep(NA, length(styr:get_start_year(data_4_model)) - 1), FIMS::get_data(data_4_model) |> - dplyr::filter(type == "index", name == "survey1") |> + dplyr::filter(type == "index", fleet == "survey1") |> dplyr::pull(value) ), fims.expected = c(rep(NA, length(styr:get_start_year(data_4_model)) - 1), report$index_expected[[3]]), @@ -395,7 +395,7 @@ fleet1_landings_results <- data.frame( observed = c( rep(NA, length(styr:get_start_year(data_4_model)) - 1), FIMS::get_data(data_4_model) |> - dplyr::filter(type == "landings", name == "fleet1") |> + dplyr::filter(type == "landings", fleet == "fleet1") |> dplyr::pull(value) ), fims.expected = c(rep(NA, length(styr:get_start_year(data_4_model)) - 1), report$landings_expected[[1]]), @@ -406,7 +406,7 @@ fleet2_landings_results <- data.frame( observed = c( rep(NA, length(styr:get_start_year(data_4_model)) - 1), FIMS::get_data(data_4_model) |> - dplyr::filter(type == "landings", name == "fleet2") |> + dplyr::filter(type == "landings", fleet == "fleet2") |> dplyr::pull(value) ), fims.expected = c(rep(NA, length(styr:get_start_year(data_4_model)) - 1), report$landings_expected[[2]]), diff --git a/content/SWFSC-sardine.qmd b/content/SWFSC-sardine.qmd index 3168597..0fe87fc 100644 --- a/content/SWFSC-sardine.qmd +++ b/content/SWFSC-sardine.qmd @@ -75,7 +75,7 @@ catch <- data.frame(year = 2005:2023, catch = c(29188.50, 53107.00, 69929.40, # geom_line() + scale_y_continuous(label = comma) -fimscatch <- tibble::tibble(type = "landings", name = "fleet1", +fimscatch <- tibble::tibble(type = "landings", fleet = "fleet1", age = NA, timing = catch$year, value = catch$catch, unit = "mt", uncertainty = 0.05) @@ -89,16 +89,20 @@ cpue <- data.frame(year = 2005:2023, obs = c(649619.0, 899635.0, 956354.0, 86328 # ggplot(cpue, aes(x = year, y = obs)) + geom_point() + geom_line() + # scale_y_continuous(label = comma) -fimsindex <- tibble::tibble(type = "index", name = "survey1", +fimsindex <- tibble::tibble(type = "index", fleet = "survey1", age = NA, timing = cpue$year, value = cpue$obs, unit = "mt", uncertainty = .3) fimsindex$unit <- "" #-----Age compositions acomps <- utils::read.csv(file.path(data_directory, "sardine_acomps.csv")) |> - dplyr::mutate(value_prop = value / Nsamp) # convert age-comp data to proportions, to match fims-demo.Rmd + dplyr::group_by(name, Yr) |> + # convert value column to proportions, to match fims-demo.Rmd + dplyr::mutate(value_prop = value / sum(value)) |> + dplyr::ungroup() -fimsage <- tibble::tibble(type = "age_comp", name = acomps$name, + +fimsage <- tibble::tibble(type = "age_comp", fleet = acomps$name, age = acomps$age, timing = acomps$Yr, value = acomps$value_prop, unit = "proportion", uncertainty = acomps$Nsamp) #fimsage$uncertainty <- 50 Leave as empirical values @@ -114,11 +118,11 @@ wtatage <- r4ss::SS_readwtatage(file.path(data_directory, "sardine_wtatage.ss_ne dplyr::filter(fleet %in% c(1, 2)) |> dplyr::mutate(fleet = ifelse(fleet == 1, "fleet1", "survey1")) |> tidyr::pivot_longer(cols = `0`:`10`, names_to = "age", values_to = "value") |> - dplyr::filter(year != 2024, !(age %in% c("9", "10"))) |> # Trim ages 9 and 10 to match age-comps + dplyr::filter(!(age %in% c("9", "10"))) |> # Trim ages 9 and 10 to match age-comps dplyr::mutate(value = value / 1000) # WAA converted from kg to mt fimswaa <- tibble::tibble( - type = "weight-at-age", - name = wtatage$fleet, + type = "weight_at_age", + fleet = wtatage$fleet, age = wtatage$age, timing = wtatage$year, value = wtatage$value, @@ -153,9 +157,9 @@ parameters <- FIMS::create_default_configurations(data = data_4_model) |> label = rep(c("inflection_point", "slope"), 2), value = c(1, 5, 1.2, 2), estimation_type = "constant", - fleet_name = rep(c("fleet1", "survey1"), each = 2) + fleet = rep(c("fleet1", "survey1"), each = 2) ), - by = c("module_name", "label", "fleet_name") + by = c("module_name", "label", "fleet") ) |> dplyr::rows_update( tibble::tibble( @@ -192,21 +196,21 @@ parameters <- FIMS::create_default_configurations(data = data_4_model) |> dplyr::rows_update( tibble::tibble( module_name = "Fleet", - fleet_name = "survey1", + fleet = "survey1", label = "log_q", value = 0, estimation_type = "constant" ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( module_name = "Fleet", - fleet_name = "fleet1", + fleet = "fleet1", label = "log_Fmort", value = log(0.2) ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( @@ -229,22 +233,22 @@ parameters <- FIMS::create_default_configurations(data = data_4_model) |> dplyr::rows_update( tibble::tibble( module_name = "Data", - fleet_name = "fleet1", + fleet = "fleet1", label = "log_sd", value = log(sqrt(log(0.01^2 + 1))), estimation_type = "fixed_effects" ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( module_name = "Data", - fleet_name = "survey1", + fleet = "survey1", label = "log_sd", value = log(sqrt(log(0.1^2 + 1))), estimation_type = "fixed_effects" ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) ``` @@ -265,7 +269,7 @@ output <- FIMS::get_estimates(fit) |> ) # Get information about the model and print a few characters to the screen -log <- FIMS::get_log_module("information") +# log <- FIMS::get_log_module("information") ``` ## Add your comparison figures @@ -411,27 +415,28 @@ ggplot2::ggsave("figures/SWFSC-sardine-recruitment.png", width = 6.8, height = 5 #Are fixed but plot for comparison's sake ##Fishery -sel_fishery <- logistic(ages, - slope = dplyr::filter(FIMS::get_estimates(fit), module_name == "Selectivity", label == "slope", module_id == 1) |> - dplyr::pull(estimated), - inflection_point = dplyr::filter(FIMS::get_estimates(fit), module_name == "Selectivity", label == "inflection_point", module_id == 1) |> - dplyr::pull(estimated)) - -names(sel_fishery) <- c("age", "fims") - -sel_fishery$ss3 <- ssres$ageselex |> - dplyr::filter(Yr == 2005, Factor == "Asel", Fleet == 1) |> - dplyr::select(as.character(0:8)) |> - t() -sel_fishery <- sel_fishery |> - reshape2::melt(id.var = "age") - -ggplot2::ggplot( - sel_fishery, - ggplot2::aes(x = age, y = value, group = variable, color = variable) -) + - ggplot2::geom_point() + - ggplot2::geom_line() +### commented out due to error "! object 'value' not found" +# sel_fishery <- logistic(ages, +# slope = dplyr::filter(FIMS::get_estimates(fit), module_name == "Selectivity", label == "slope", module_id == 1) |> +# dplyr::pull(estimated), +# inflection_point = dplyr::filter(FIMS::get_estimates(fit), module_name == "Selectivity", label == "inflection_point", module_id == 1) |> +# dplyr::pull(estimated)) + +# names(sel_fishery) <- c("age", "fims") + +# sel_fishery$ss3 <- ssres$ageselex |> +# dplyr::filter(Yr == 2005, Factor == "Asel", Fleet == 1) |> +# dplyr::select(as.character(0:8)) |> +# t() +# sel_fishery <- sel_fishery |> +# reshape2::melt(id.var = "age") + +# ggplot2::ggplot( +# sel_fishery, +# ggplot2::aes(x = age, y = value, group = variable, color = variable) +# ) + +# ggplot2::geom_point() + +# ggplot2::geom_line() #-----------Survey sel_survey <- logistic(ages, diff --git a/content/advanced-features.qmd b/content/advanced-features.qmd index ccf89bf..a812236 100644 --- a/content/advanced-features.qmd +++ b/content/advanced-features.qmd @@ -83,14 +83,14 @@ estimate_recdevs <- TRUE # Set up a FIMS model without wrapper functions -fishery_catch <- FIMS::m_landings(data_4_model, "fleet1") -fishery_agecomp <- FIMS::m_agecomp(data_4_model, "fleet1") -survey_index2 <- FIMS::m_index(data_4_model, "survey2") -survey_agecomp2 <- FIMS::m_agecomp(data_4_model, "survey2") -survey_index3 <- FIMS::m_index(data_4_model, "survey3") -survey_agecomp3 <- FIMS::m_agecomp(data_4_model, "survey3") -survey_index6 <- FIMS::m_index(data_4_model, "survey6") -survey_agecomp6 <- FIMS::m_agecomp(data_4_model, "survey6") +fishery_catch <- FIMS::model_landings(data_4_model, "fleet1") +fishery_agecomp <- FIMS::model_age_comp(data_4_model, "fleet1") +survey_index2 <- FIMS::model_index(data_4_model, "survey2") +survey_agecomp2 <- FIMS::model_age_comp(data_4_model, "survey2") +survey_index3 <- FIMS::model_index(data_4_model, "survey3") +survey_agecomp3 <- FIMS::model_age_comp(data_4_model, "survey3") +survey_index6 <- FIMS::model_index(data_4_model, "survey6") +survey_agecomp6 <- FIMS::model_age_comp(data_4_model, "survey6") # need to think about how to deal with multiple fleets - only using 1 fleet for now # TODO: FIMS now supports multiple fishing fleets. # We can test this feature using the case study to evaluate its functionality. @@ -108,7 +108,7 @@ purrr::walk( dplyr::filter( .data = get_data(data_4_model), type == "age_comp", - name %in% "fleet1" + fleet %in% "fleet1" ) |> dplyr::pull(uncertainty) )[x] @@ -139,7 +139,7 @@ for (y in seq(n_years)) { # Log-transform OM fishing mortality fish_fleet$log_Fmort[y]$value <- log(pkfitfinal$rep$F[y]) } -fish_fleet$log_Fmort$set_all_estimable(TRUE) +fish_fleet$log_Fmort$set_estimation_types(c("fixed_effects")) fish_fleet$log_q[1]$value <- log(1.0) # why is this length two in Chris' case study? fish_fleet$log_q[1]$estimation_type$set("constant") @@ -156,7 +156,7 @@ for (y in seq(n_years)) { # Compute lognormal SD from OM coefficient of variation (CV) fish_fleet_index_distribution$log_sd[y]$value <- log(fimsdat$cattot_log_sd[y]) } -fish_fleet_index_distribution$log_sd$set_all_estimable(FALSE) +fish_fleet_index_distribution$log_sd$set_estimation_types(c("constant")) # Set Data using the IDs from the modules defined above fish_fleet_index_distribution$set_observed_data(fish_fleet$GetObservedLandingsDataID()) fish_fleet_index_distribution$set_distribution_links("data", fish_fleet$log_index_expected$get_id()) @@ -181,7 +181,7 @@ purrr::walk( dplyr::filter( .data = get_data(data_4_model), type == "age_comp", - name %in% "survey2" + fleet %in% "survey2" ) |> dplyr::pull(uncertainty) )[x] @@ -209,7 +209,7 @@ for (y in seq_along(years)) { # Set very low survey fishing mortality survey2_fleet$log_Fmort[y]$value <- -200 } -survey2_fleet$log_Fmort$set_all_estimable(FALSE) +survey2_fleet$log_Fmort$set_estimation_types(c("constant")) survey2_fleet$log_q[1]$value <- parfinal$log_q2_mean survey2_fleet$log_q[1]$estimation_type$set("fixed_effects") survey2_fleet$SetSelectivityID(survey2_selex$get_id()) @@ -232,7 +232,7 @@ for (y in which(temporary[, "index"] != -999)) { ) } rm(temporary) -survey2_fleet_index_distribution$log_sd$set_all_estimable(FALSE) +survey2_fleet_index_distribution$log_sd$set_estimation_types(c("constant")) # Set Data using the IDs from the modules defined above survey2_fleet_index_distribution$set_observed_data(survey2_fleet$GetObservedIndexDataID()) survey2_fleet_index_distribution$set_distribution_links("data", survey2_fleet$log_index_expected$get_id()) @@ -256,7 +256,7 @@ purrr::walk( dplyr::filter( .data = get_data(data_4_model), type == "age_comp", - name %in% "survey3" + fleet %in% "survey3" ) |> dplyr::pull(uncertainty) )[x] @@ -278,7 +278,7 @@ for (y in seq_along(years)) { # Set very low survey fishing mortality survey3_fleet$log_Fmort[y]$value <- -200 } -survey3_fleet$log_Fmort$set_all_estimable(FALSE) +survey3_fleet$log_Fmort$set_estimation_types(c("constant")) survey3_fleet$log_q[1]$value <- parfinal$log_q3_mean survey3_fleet$log_q[1]$estimation_type$set("fixed_effects") survey3_fleet$SetSelectivityID(survey3_selex$get_id()) @@ -302,7 +302,7 @@ for (y in which(temporary[, "index"] != -999)) { ) } rm(temporary) -survey3_fleet_index_distribution$log_sd$set_all_estimable(FALSE) +survey3_fleet_index_distribution$log_sd$set_estimation_types(c("constant")) # Set Data using the IDs from the modules defined above survey3_fleet_index_distribution$set_observed_data(survey3_fleet$GetObservedIndexDataID()) survey3_fleet_index_distribution$set_distribution_links("data", survey3_fleet$log_index_expected$get_id()) @@ -326,7 +326,7 @@ purrr::walk( dplyr::filter( .data = get_data(data_4_model), type == "age_comp", - name %in% "survey6" + fleet %in% "survey6" ) |> dplyr::pull(uncertainty) )[x] @@ -354,7 +354,7 @@ for (y in seq_along(years)) { # Set very low survey fishing mortality survey6_fleet$log_Fmort[y]$value <- -200 } -survey6_fleet$log_Fmort$set_all_estimable(FALSE) +survey6_fleet$log_Fmort$set_estimation_types(c("constant")) survey6_fleet$log_q[1]$value <- parfinal$log_q6 survey6_fleet$log_q[1]$estimation_type$set("fixed_effects") survey6_fleet$SetSelectivityID(survey6_selex$get_id()) @@ -377,7 +377,7 @@ for (y in which(temporary[, "index"] != -999)) { ) } rm(temporary) -survey6_fleet_index_distribution$log_sd$set_all_estimable(FALSE) +survey6_fleet_index_distribution$log_sd$set_estimation_types(c("constant")) # Set Data using the IDs from the modules defined above survey6_fleet_index_distribution$set_observed_data(survey6_fleet$GetObservedIndexDataID()) survey6_fleet_index_distribution$set_distribution_links("data", survey6_fleet$log_index_expected$get_id()) @@ -403,7 +403,7 @@ recruitment$log_devs$resize(n_years - 1) for (y in seq(n_years - 1)) { recruitment$log_devs[y]$value <- parfinal$dev_log_recruit[y + 1] } -recruitment$log_devs$set_all_random(TRUE) +recruitment$log_devs$set_estimation_types(c("random_effects")) recruitment_distribution <- methods::new(DnormDistribution) # set up logR_sd using the normal log_sd parameter # logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is @@ -415,10 +415,10 @@ recruitment_distribution$log_sd[1]$value <- log(parfinal$sigmaR) # TODO: should be estimated b/c it is a random effect # TODO: sigma_R doesn't have a variable map yet so you cannot set a prior recruitment_distribution$log_sd[1]$estimation_type$set("constant") -recruitment_distribution$x$resize(n_years - 1) +recruitment_distribution$observed_values$resize(n_years - 1) recruitment_distribution$expected_values$resize(n_years - 1) for (i in seq(n_years - 1)) { - recruitment_distribution$x[i]$value <- 0 + recruitment_distribution$observed_values[i]$value <- 0 recruitment_distribution$expected_values[i]$value <- 0 } recruitment_distribution$set_distribution_links("random_effects", recruitment$log_devs$get_id()) @@ -428,7 +428,7 @@ recruitment_distribution$set_distribution_links("random_effects", recruitment$lo waa <- pkinput$dat$wt_srv1 waa <- rbind(waa, waa[1, ]) ewaa_growth <- methods::new(EWAAGrowth) -ewaa_growth$n_years$set(get_n_years(data_4_model) + 1) +ewaa_growth$n_years$set(get_n_years(data_4_model)) ewaa_growth$ages$resize(n_ages) purrr::walk( seq_along(ages), @@ -460,13 +460,13 @@ population$log_M$resize(n_years * n_ages) for (i in seq(n_years * n_ages)) { population$log_M[i]$value <- tmpM[i] } -population$log_M$set_all_estimable(FALSE) +population$log_M$set_estimation_types(c("constant")) population$log_init_naa$resize(n_ages) initNAA <- c(log(pkfitfinal$rep$recruit[1]), log(pkfitfinal$rep$initN)) + log(1e9) for (i in seq(n_ages)) { population$log_init_naa[i]$value <- initNAA[i] } -population$log_init_naa$set_all_estimable(FALSE)# NOTE: fixing at ASAP estimates to test SSB calculations +population$log_init_naa$set_estimation_types(c("constant"))# NOTE: fixing at ASAP estimates to test SSB calculations population$n_ages$set(n_ages) population$ages$resize(n_ages) purrr::walk( @@ -507,7 +507,7 @@ opt <- with(obj, nlminb( ) )) FIMS::set_fixed(opt$par) -fims_finalized <- caa$get_output(do_sd_report = FALSE) +fims_finalized <- caa$get_output() max(abs(obj$gr())) # from Cole, can use TMBhelper::fit_tmb to get val to <1e-10 # FIMS after estimation rep2 <- obj$report(obj$env$last.par.best) @@ -744,14 +744,14 @@ data_4_model <- prepare_pollock_data( ages = ages, n_ages = n_ages ) -fishery_catch <- FIMS::m_landings(data_4_model, "fleet1") -fishery_agecomp <- FIMS::m_agecomp(data_4_model, "fleet1") -survey_index2 <- FIMS::m_index(data_4_model, "survey2") -survey_agecomp2 <- FIMS::m_agecomp(data_4_model, "survey2") -survey_index3 <- FIMS::m_index(data_4_model, "survey3") -survey_agecomp3 <- FIMS::m_agecomp(data_4_model, "survey3") -survey_index6 <- FIMS::m_index(data_4_model, "survey6") -survey_agecomp6 <- FIMS::m_agecomp(data_4_model, "survey6") +fishery_catch <- FIMS::model_landings(data_4_model, "fleet1") +fishery_agecomp <- FIMS::model_age_comp(data_4_model, "fleet1") +survey_index2 <- FIMS::model_index(data_4_model, "survey2") +survey_agecomp2 <- FIMS::model_age_comp(data_4_model, "survey2") +survey_index3 <- FIMS::model_index(data_4_model, "survey3") +survey_agecomp3 <- FIMS::model_age_comp(data_4_model, "survey3") +survey_index6 <- FIMS::model_index(data_4_model, "survey6") +survey_agecomp6 <- FIMS::model_age_comp(data_4_model, "survey6") # need to think about how to deal with multiple fleets - only using 1 fleet for now # TODO: FIMS now supports multiple fishing fleets. # We can test this feature using the case study to evaluate its functionality. @@ -769,7 +769,7 @@ purrr::walk( dplyr::filter( .data = get_data(data_4_model), type == "age_comp", - name %in% "fleet1" + fleet %in% "fleet1" ) |> dplyr::pull(value) )[x] @@ -792,15 +792,15 @@ fish_selex$slope_desc[1]$estimation_type$set(estimate_fish_selex) ## create fleet object fish_fleet <- methods::new(Fleet) -fish_fleet$nages$set(n_ages) -fish_fleet$nyears$set(n_years) +fish_fleet$n_ages$set(n_ages) +fish_fleet$n_years$set(n_years) fish_fleet$log_Fmort$resize(n_years) for (y in seq(n_years)) { # Log-transform OM fishing mortality fish_fleet$log_Fmort[y]$value <- log(pkfitfinal$rep$F[y]) } -fish_fleet$log_Fmort$set_all_estimable(TRUE) +fish_fleet$log_Fmort$set_estimation_types(c("fixed_effects")) fish_fleet$log_q[1]$value <- 0 # why is this length two in Chris' case study? fish_fleet$log_q[1]$estimation_type$set("constant") @@ -815,9 +815,9 @@ fish_fleet_index_distribution <- methods::new(DlnormDistribution) fish_fleet_index_distribution$log_sd$resize(n_years) for (y in seq(n_years)) { # Compute lognormal SD from OM coefficient of variation (CV) - fish_fleet_index_distribution$log_sd[y]$value <- log(landings$uncertainty[y]) + fish_fleet_index_distribution$log_sd[y]$value <- log(fimsdat$cattot_log_sd[y]) } -fish_fleet_index_distribution$log_sd$set_all_estimable(FALSE) +fish_fleet_index_distribution$log_sd$set_estimation_types("constant") # Set Data using the IDs from the modules defined above fish_fleet_index_distribution$set_observed_data(fish_fleet$GetObservedIndexDataID()) fish_fleet_index_distribution$set_distribution_links("data", fish_fleet$log_index_expected$get_id()) @@ -842,7 +842,7 @@ purrr::walk( dplyr::filter( .data = get_data(data_4_model), type == "age_comp", - name %in% "survey2" + fleet %in% "survey2" ) |> dplyr::pull(value) )[x] @@ -863,8 +863,8 @@ survey2_selex$slope_desc[1]$value <- exp(parfinal$log_slp2_srv2) survey2_selex$slope_desc[1]$estimation_type$set("constant") survey2_fleet <- methods::new(Fleet) -survey2_fleet$nages$set(n_ages) -survey2_fleet$nyears$set(n_years) +survey2_fleet$n_ages$set(n_ages) +survey2_fleet$n_years$set(n_years) survey2_fleet$log_q[1]$value <- parfinal$log_q2_mean survey2_fleet$log_q[1]$estimation_type$set("fixed_effects") survey2_fleet$SetSelectivityID(survey2_selex$get_id()) @@ -876,9 +876,9 @@ survey2_fleet_index_distribution <- methods::new(DlnormDistribution) survey2_fleet_index_distribution$log_sd$resize(n_years) for (y in seq(n_years)) { # Compute lognormal SD from OM coefficient of variation (CV) - survey2_fleet_index_distribution$log_sd[y]$value <- log(index2$uncertainty)[y] + survey2_fleet_index_distribution$log_sd[y]$value <- log(dplyr::filter(get_data(data_4_model), fleet == "survey2", type == "index")$uncertainty)[y] } -survey2_fleet_index_distribution$log_sd$set_all_estimable(FALSE) +survey2_fleet_index_distribution$log_sd$set_estimation_types("constant") # Set Data using the IDs from the modules defined above survey2_fleet_index_distribution$set_observed_data(survey2_fleet$GetObservedIndexDataID()) survey2_fleet_index_distribution$set_distribution_links("data", survey2_fleet$log_index_expected$get_id()) @@ -902,7 +902,7 @@ purrr::walk( dplyr::filter( .data = get_data(data_4_model), type == "age_comp", - name %in% "survey3" + fleet %in% "survey3" ) |> dplyr::pull(value) )[x] @@ -917,8 +917,8 @@ survey3_selex$slope[1]$value <- exp(parfinal$log_slp1_srv3) survey3_selex$slope[1]$estimation_type$set(estimate_survey_selex) survey3_fleet <- methods::new(Fleet) -survey3_fleet$nages$set(n_ages) -survey3_fleet$nyears$set(n_years) +survey3_fleet$n_ages$set(n_ages) +survey3_fleet$n_years$set(n_years) survey3_fleet$log_q[1]$value <- parfinal$log_q3_mean survey3_fleet$log_q[1]$estimation_type$set("fixed_effects") survey3_fleet$SetSelectivityID(survey3_selex$get_id()) @@ -931,9 +931,9 @@ survey3_fleet_index_distribution <- methods::new(DlnormDistribution) survey3_fleet_index_distribution$log_sd$resize(n_years) for (y in seq(n_years)) { # Compute lognormal SD from OM coefficient of variation (CV) - survey3_fleet_index_distribution$log_sd[y]$value <- log(index3$uncertainty)[y] + survey3_fleet_index_distribution$log_sd[y]$value <- log(1) } -survey3_fleet_index_distribution$log_sd$set_all_estimable(FALSE) +survey3_fleet_index_distribution$log_sd$set_estimation_types("constant") # Set Data using the IDs from the modules defined above survey3_fleet_index_distribution$set_observed_data(survey3_fleet$GetObservedIndexDataID()) survey3_fleet_index_distribution$set_distribution_links("data", survey3_fleet$log_index_expected$get_id()) @@ -957,7 +957,7 @@ purrr::walk( dplyr::filter( .data = get_data(data_4_model), type == "age_comp", - name %in% "survey6" + fleet %in% "survey6" ) |> dplyr::pull(value) )[x] @@ -978,8 +978,8 @@ survey6_selex$slope_desc[1]$value <- exp(parfinal$log_slp2_srv6) survey6_selex$slope_desc[1]$estimation_type$set(estimate_survey_selex) survey6_fleet <- methods::new(Fleet) -survey6_fleet$nages$set(n_ages) -survey6_fleet$nyears$set(n_years) +survey6_fleet$n_ages$set(n_ages) +survey6_fleet$n_years$set(n_years) survey6_fleet$log_q[1]$value <- parfinal$log_q6 survey6_fleet$log_q[1]$estimation_type$set("fixed_effects") survey6_fleet$SetSelectivityID(survey6_selex$get_id()) @@ -991,9 +991,9 @@ survey6_fleet_index_distribution <- methods::new(DlnormDistribution) survey6_fleet_index_distribution$log_sd$resize(n_years) for (y in seq(n_years)) { # Compute lognormal SD from OM coefficient of variation (CV) - survey6_fleet_index_distribution$log_sd[y]$value <- log(index6$uncertainty)[y] + survey6_fleet_index_distribution$log_sd[y]$value <- log(1)[y] } -survey6_fleet_index_distribution$log_sd$set_all_estimable(FALSE) +survey6_fleet_index_distribution$log_sd$set_estimation_types("constant") # Set Data using the IDs from the modules defined above survey6_fleet_index_distribution$set_observed_data(survey6_fleet$GetObservedIndexDataID()) survey6_fleet_index_distribution$set_distribution_links("data", survey6_fleet$log_index_expected$get_id()) @@ -1015,24 +1015,23 @@ recruitment$log_rzero[1]$estimation_type$set("fixed_effects") recruitment$logit_steep[1]$value <- -log(1.0 - .99999) + log(.99999 - 0.2) recruitment$logit_steep[1]$estimation_type$set("constant") -recruitment$nyears$set(n_years - 1) +recruitment$n_years$set(n_years - 1) recruitment$log_devs$resize(n_years - 1) for (y in seq(n_years - 1)) { recruitment$log_devs[y]$value <- parfinal$dev_log_recruit[y+1] } -recruitment$log_devs$set_all_estimable(estimate_recdevs) -recruitment$log_devs$set_all_random(TRUE) +recruitment$log_devs$set_estimation_types("random_effects") recruitment_distribution <- methods::new(DnormDistribution) # set up logR_sd using the normal log_sd parameter # logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is # taken before the likelihood calculation -recruitment_distribution$log_sd <- methods::new(ParameterVector, 1) +recruitment_distribution$log_sd <- methods::new(VariableVector, 1) recruitment_distribution$log_sd[1]$value <- log(parfinal$sigmaR) recruitment_distribution$log_sd[1]$estimation_type$set("constant") -recruitment_distribution$x$resize(n_years - 1) +recruitment_distribution$observed_values$resize(n_years - 1) recruitment_distribution$expected_values$resize(n_years - 1) for (i in seq(n_years - 1)) { - recruitment_distribution$x[i]$value <- 0 + recruitment_distribution$observed_values[i]$value <- 0 recruitment_distribution$expected_values[i]$value <- 0 } recruitment_distribution$set_distribution_links("random_effects", recruitment$log_devs$get_id()) @@ -1042,7 +1041,7 @@ recruitment_distribution$set_distribution_links("random_effects", recruitment$lo waa <- pkinput$dat$wt_srv1 waa <- rbind(waa, waa[1, ]) ewaa_growth <- methods::new(EWAAGrowth) -ewaa_growth$n_years$set(get_n_years(data_4_model) + 1) +ewaa_growth$n_years$set(get_n_years(data_4_model)) ewaa_growth$ages$resize(n_ages) purrr::walk( seq_along(ages), @@ -1074,22 +1073,21 @@ population$log_M$resize(n_years * n_ages) for (i in seq(n_years * n_ages)) { population$log_M[i]$value <- tmpM[i] } -population$log_M$set_all_estimable(FALSE) +population$log_M$set_estimation_types("constant") population$log_init_naa$resize(n_ages) initNAA <- c(log(pkfitfinal$rep$recruit[1]), log(pkfitfinal$rep$initN)) + log(1e9) for (i in seq(n_ages)) { population$log_init_naa[i]$value <- initNAA[i] } -population$log_init_naa$set_all_estimable(FALSE)# NOTE: fixing at ASAP estimates to test SSB calculations -population$nages$set(n_ages) +population$log_init_naa$set_estimation_types("constant")# NOTE: fixing at ASAP estimates to test SSB calculations +population$n_ages$set(n_ages) population$ages$resize(n_ages) purrr::walk( seq_along(ages), \(x) population$ages$set(x - 1, ages[x]) ) -population$nfleets$set(4) -population$nyears$set(n_years) -population$nseasons$set(1) +population$n_fleets$set(4) +population$n_years$set(n_years) population$SetMaturityID(maturity$get_id()) population$SetGrowthID(ewaa_growth$get_id()) population$SetRecruitmentID(recruitment$get_id()) diff --git a/content/data_files/opaka_model.RDS b/content/data_files/opaka_model.RDS index 6aba344..800d4aa 100644 Binary files a/content/data_files/opaka_model.RDS and b/content/data_files/opaka_model.RDS differ diff --git a/content/figures/SEFSC_scamp_caa_fleet1.png b/content/figures/SEFSC_scamp_caa_fleet1.png index 080989c..ff180ed 100644 Binary files a/content/figures/SEFSC_scamp_caa_fleet1.png and b/content/figures/SEFSC_scamp_caa_fleet1.png differ diff --git a/content/figures/SEFSC_scamp_caa_fleet2.png b/content/figures/SEFSC_scamp_caa_fleet2.png index 9fab8b7..12f7621 100644 Binary files a/content/figures/SEFSC_scamp_caa_fleet2.png and b/content/figures/SEFSC_scamp_caa_fleet2.png differ diff --git a/content/figures/SEFSC_scamp_caa_survey.png b/content/figures/SEFSC_scamp_caa_survey.png index 15fd5b8..a5192ce 100644 Binary files a/content/figures/SEFSC_scamp_caa_survey.png and b/content/figures/SEFSC_scamp_caa_survey.png differ diff --git a/content/figures/SEFSC_scamp_selex.png b/content/figures/SEFSC_scamp_selex.png index d33461f..93497bf 100644 Binary files a/content/figures/SEFSC_scamp_selex.png and b/content/figures/SEFSC_scamp_selex.png differ diff --git a/content/figures/SEFSC_scamp_tseries_F.png b/content/figures/SEFSC_scamp_tseries_F.png index 717013d..9cd7797 100644 Binary files a/content/figures/SEFSC_scamp_tseries_F.png and b/content/figures/SEFSC_scamp_tseries_F.png differ diff --git a/content/figures/SEFSC_scamp_tseries_fits.png b/content/figures/SEFSC_scamp_tseries_fits.png index 5a1eee8..cdd33e9 100644 Binary files a/content/figures/SEFSC_scamp_tseries_fits.png and b/content/figures/SEFSC_scamp_tseries_fits.png differ diff --git a/content/figures/SEFSC_scamp_tseries_popn.png b/content/figures/SEFSC_scamp_tseries_popn.png index 9b3469d..95ae13c 100644 Binary files a/content/figures/SEFSC_scamp_tseries_popn.png and b/content/figures/SEFSC_scamp_tseries_popn.png differ diff --git a/content/pacific-hake.qmd b/content/pacific-hake.qmd index e29d55e..5097578 100644 --- a/content/pacific-hake.qmd +++ b/content/pacific-hake.qmd @@ -84,7 +84,7 @@ stock_synthesis_data <- r4ss::ss3_data_to_fims( ) |> dplyr::mutate( uncertainty = ifelse( - test = type == "age" & name == "Acoustic_Survey", + test = type == "age" & fleet == "Acoustic_Survey", uncertainty * 30, uncertainty ) @@ -133,7 +133,7 @@ data_4_model # Summary of the data types available FIMS::get_data(data_4_model) |> dplyr::filter(value != -999) |> - dplyr::group_by(type, name) |> + dplyr::group_by(type, fleet) |> dplyr::count() ``` @@ -283,10 +283,10 @@ parameters <- FIMS::create_default_configurations(data = data_4_model) |> dplyr::rows_update( tibble::tibble( module_name = "Selectivity", - fleet_name = c("Fishery", "Acoustic_Survey"), + fleet = c("Fishery", "Acoustic_Survey"), module_type = c(selectivity_used, "Logistic") ), - by = c("module_name", "fleet_name") + by = c("module_name", "fleet") ) |> FIMS::create_default_parameters(data = data_4_model) |> tidyr::unnest(data) |> @@ -325,45 +325,45 @@ parameters <- FIMS::create_default_configurations(data = data_4_model) |> dplyr::rows_update( tibble::tibble( module_name = "Selectivity", - fleet_name = rep(c("Fishery", "Acoustic_Survey"), each = 2), + fleet = rep(c("Fishery", "Acoustic_Survey"), each = 2), label = rep(c("inflection_point", "slope"), 2), # TODO: play with these start values value = c(2.05, 3.17, 4.4, 3.3), estimation_type = "fixed_effects" ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) |> dplyr::rows_update( tibble::tibble( module_name = "Fleet", - fleet_name = "Acoustic_Survey", + fleet = "Acoustic_Survey", label = c("log_q"), value = log(0.832), estimation_type = "fixed_effects" ), - by = c("module_name", "fleet_name", "label") + by = c("module_name", "fleet", "label") ) -fleet_names <- get_fleets(data_4_model) +fleets <- get_fleets(data_4_model) # Initialize lists to store fleet-related objects fleet <- fleet_selectivity <- fleet_landings <- fleet_landings_distribution <- fleet_index <- fleet_index_distribution <- fleet_age_comp <- fleet_agecomp_distribution <- - vector("list", length(fleet_names)) + vector("list", length(fleets)) -for (i in seq_along(fleet_names)) { +for (i in seq_along(fleets)) { # Selectivity fleet_selectivity[[i]] <- FIMS:::initialize_selectivity( parameters = parameters, data = data_4_model, - fleet_name = fleet_names[i] + fleet = fleets[i] ) # Time-varying selectivity for the fishery if (i == 1 & time_varying_selectivity) { selectivity_type <- parameters |> - dplyr::filter(module_name == "Selectivity", fleet_name == "Fishery") |> + dplyr::filter(module_name == "Selectivity", fleet == "Fishery") |> dplyr::slice(1) |> dplyr::pull(module_type) if (selectivity_type == "Logistic"){ @@ -397,11 +397,11 @@ for (i in seq_along(fleet_names)) { selectivity = fleet_selectivity[[i]]$get_id() ) fleet_types <- get_data(data_4_model) |> - dplyr::filter(name == fleet_names[i]) |> + dplyr::filter(fleet == fleets[i]) |> dplyr::pull(type) |> unique() data_distribution_names_for_fleet_i <- parameters |> - dplyr::filter(fleet_name == fleet_names[i] & distribution_type == "Data") |> + dplyr::filter(fleet == fleets[i] & distribution_type == "Data") |> dplyr::pull(module_type) # Landings @@ -410,7 +410,7 @@ for (i in seq_along(fleet_names)) { # Initialize landings module for the current fleet fleet_landings[[i]] <- FIMS:::initialize_landings( data = data_4_model, - fleet_name = fleet_names[i] + fleet = fleets[i] ) # Add the module ID for the initialized landings to the list of fleet module IDs fleet_module_ids <- c( @@ -425,7 +425,7 @@ for (i in seq_along(fleet_names)) { # Initialize index module for the current fleet fleet_index[[i]] <- FIMS:::initialize_index( data = data_4_model, - fleet_name = fleet_names[i] + fleet = fleets[i] ) fleet_module_ids <- c( fleet_module_ids, @@ -439,7 +439,7 @@ for (i in seq_along(fleet_names)) { # Initialize age composition module for the current fleet fleet_age_comp[[i]] <- FIMS:::initialize_comp( data = data_4_model, - fleet_name = fleet_names[i], + fleet = fleets[i], type = "AgeComp" ) fleet_module_ids <- c( @@ -451,13 +451,13 @@ for (i in seq_along(fleet_names)) { fleet[[i]] <- FIMS:::initialize_fleet( parameters = parameters, data = data_4_model, - fleet_name = fleet_names[i], + fleet = fleets[i], linked_ids = fleet_module_ids ) # Fleet uncertainty fleet_sd_input <- parameters |> - dplyr::filter(fleet_name == fleet_names[i] & label == "log_sd") |> + dplyr::filter(fleet == fleets[i] & label == "log_sd") |> dplyr::mutate( label = "sd", value = exp(value) @@ -511,7 +511,7 @@ recruitment$log_devs$resize(get_n_years(data_4_model)) for (y in seq(recruitment$log_devs$size())) { recruitment$log_devs[y]$value <- 0 } -recruitment$log_devs$set_all_estimable(TRUE) +recruitment$log_devs$set_estimation_types(c("constant")) recruitment$n_years$set(get_n_years(data_4_model)) recruitment_distribution <- methods::new(DnormDistribution) # TODO: check with Andrea about log space here @@ -532,7 +532,8 @@ recruitment_distribution$set_distribution_links( "random_effects", recruitment$log_devs$get_id() ) -recruitment$log_devs$set_all_random(TRUE) +recruitment$log_devs$set_estimation_types(c("random_effects")) + # recruitment$log_devs[46]$estimation_type$set("constant") # recruitment$log_devs[47]$estimation_type$set("constant") # recruitment$log_devs[48]$estimation_type$set("constant") @@ -790,7 +791,7 @@ posteriors_labeled <- tidyr::pivot_longer( } else { dplyr::row_number(label) }, - fleet_name = dplyr::if_else(fleet_number == 1, "Fishery", "Acoustic_Survey"), + fleet = dplyr::if_else(fleet_number == 1, "Fishery", "Acoustic_Survey"), words = gsub("_", " ", label) ) |> dplyr::ungroup(), @@ -835,7 +836,7 @@ ggplot2::ggplot( alpha = 0.6, .width = c(.66, .95) # Shows both 66% and 95% intervals ) + - ggplot2::facet_wrap(fleet_name ~ words, scales = "free_x") + + ggplot2::facet_wrap(fleet ~ words, scales = "free_x") + ggplot2::labs(title = "Posterior Distribution", x = "Parameter Value", y = "Density") + ggplot2::theme_minimal() ```