In Bai-Li-NOAA/IFA4EBFM: Interdisciplinary Forecasting Approach for Ecosystem-Based Fisheries Management

knitr::opts_chunk$set(eval = FALSE, echo = TRUE) # Skip the vignette when eval = FALSE

Simulation scenarios {.tabset}

Environmental drivers have been incorporated through either forcing function (producer) or egg production function. When applying forcing function, environmental drivers have been linked with phytoplankton. When applying egg production function, only AMO without lag and with lag of one year have been linked with menhaden-like species.

Environmental drivers

AMO0: Atlantic Multidecadal Oscillation without lag
AMO1: Atlantic Multidecadal Oscillation with lag of 1 year
PCP: precipitation (average value from VA, MD, CT, RI, and MA)
PDSI: Palmer drought severity index (average value from VA, MD, CT, RI, and MA)
SPLINE11: Primary producer spline 11
SPLINE33: Primary producer spline 33
SST: Kaplan sea surface temperature

{-}

devtools::load_all()

required_pkg <- c(
  "here", "ggplot2", "ggcorrplot", "RColorBrewer",
  "tidyverse", "reshape2"
)

pkg_to_install <- required_pkg[!(required_pkg %in%
                                   installed.packages()[, "Package"])]
if (length(pkg_to_install)) install.packages(pkg_to_install)

lapply(required_pkg, library, character.only = TRUE)

Biomass over time {.tabset}

ewe_year <- 1985:2017
# EwE functional groups
functional_groups = c(
  "StripedBass0",
  "StripedBass2_5",
  "StripedBass6",
  "AtlanticMenhaden0",
  "AtlanticMenhaden1",
  "AtlanticMenhaden2",
  "AtlanticMenhaden3",
  "AtlanticMenhaden4",
  "AtlanticMenhaden5",
  "AtlanticMenhaden6",
  "SpinyDogfish",
  "BluefishJuvenile",
  "BluefishAdult",
  "WeakfishJuvenile",
  "WeakfishAdult",
  "AtlanticHerring0_1",
  "AtlanticHerring2",
  "Anchovies",
  "Benthos",
  "Zooplankton",
  "Phytoplankton",
  "Detritus"
)

scenario_code <- paste0(
  "ecosim_", 
  c(
    "base_run", 
    "forcing_pcp_egg_amo1",
    "forcing_pdsi_egg_amo1",
    "forcing_spline11_egg_amo1",
    "forcing_spline33_egg_amo1",
    "forcing_sst_egg_amo1"
  ))

scenario_name <- c(
  "Base Case", 
  "PCP+AMO1",
  "PDSI+AMO1",
  "SPLINE11+AMO1",
  "SPLINE33+AMO1",
  "SST+AMO1"
)

labels <- c(
  "Striped Bass 0", "Striped Bass 2-5", "Striped Bass 6+",
  "Atlantic Menhaden 0", "Atlantic Menhaden 1",
  "Atlantic Menhaden 2", "Atlantic Menhaden 3",
  "Atlantic Menhaden 4", "Atlantic Menhaden 5",
  "Atlantic Menhaden 6+",
  "Spiny Dogfish", "Bluefish juv", "Bluefish adult",
  "Weakfish juv", "Weakfish adult",
  "Atlantic Herring 0-1", "Atlantic Herring 2+",
  "Anchovies", "Benthos", "Zooplkanton", "Phytoplankton", "Detritus"
)

reshape_ewe_data <- function(file_path, file_name, skip_nrows, functional_groups){

  for (id in 1:length(file_path)){
    ewe_temp <- read_ewe_output(
      file_path = file_path[id],
      file_names = file_name,
      skip_nrows = skip_nrows,
      plot = FALSE,
      figure_titles = NULL,
      functional_groups = functional_groups,
      figure_colors = NULL
    )

    time_id <- seq(1, nrow(ewe_temp[[1]]), by = 12)

    ewe_output <- list(ewe_temp[[1]][time_id, ])
    ewe_output[[1]]$Year <- ewe_year

    names(ewe_output[[1]])[2:ncol(ewe_output[[1]])] <- labels
    ewe_output_melt <- reshape2::melt(ewe_output[[1]], id.vars = "Year")
    ewe_output_melt$Scenario <- scenario_name[id]

    if (id==1) {
      ewe_data <- ewe_output_melt
    } else{
      ewe_data <- rbind(ewe_data, ewe_output_melt)
    }
  }

  return(ewe_data)
}

# Reshape biomass data from different scenarios
file_path <- paste0(here::here("data", "ewe", "7ages_newsim_final", "ewe7ages_ecosim_final", scenario_code))
ewe_data <- reshape_ewe_data(file_path = file_path,
                   file_name = "biomass_monthly.csv",
                   skip_nrows = 8,
                   functional_groups = functional_groups)

Biomass trends from cases that use forcing function and AMO index through egg production

subscenarios <- paste0(
  "ecosim_", 
  c(
    "base_run",
    "forcing_pcp_egg_amo1",
    "forcing_pdsi_egg_amo1",
    "forcing_spline11_egg_amo1",
    "forcing_spline33_egg_amo1",
    "forcing_sst_egg_amo1"
  ))
subdata <- ewe_data[which(ewe_data$Scenario %in% scenario_name[scenario_code %in% subscenarios]), ]
biomass <- ggplot(subdata, aes(x = Year, y = value, color = Scenario)) +
  geom_line() +
  labs(
    title = "Comparison between scenarios",
    x = "Year",
    y = "Biomass (million mt)"
  ) +
  theme_bw()

biomass + facet_wrap(~variable, scales = "free_y") + scale_color_brewer(palette="Dark2")

Total biomass of menhaden-like species from cases that use forcing function and AMO index through egg production

ss_biomass <- subdata[grep("Atlantic Menhaden", subdata$variable), ]
ss_biomass <- aggregate(value~Year+Scenario, ss_biomass, sum)

total_biomass <- ggplot(ss_biomass, aes(x = Year, y = value, color = Scenario)) +
  geom_line() +
  labs(
    title = "Comparison between scenarios",
    x = "Year",
    y = "Total Biomass (million mt)"
  ) +
  theme_bw()

total_biomass + scale_color_brewer(palette="Dark2")

file_path <- here::here("data", "ewe", "7ages_newsim_final", "ewe7ages_ecosim_final", "ecosim_forcing_pdsi_egg_amo1")
source(here::here("Rscript", "load_indices.R"))

# Load ewe value data
file_path <- paste0(here::here("data", "ewe", "7ages_newsim_final", "ewe7ages_ecosim_final", scenario_code))
ewe_catchability <- read.csv(here::here("data", "ewe", "7ages_newsim_final", "qatage.csv"))
ewe_faa <- read.csv(here::here("data", "ewe", "7ages_newsim_final", "fatage.csv"))

# Load striped bass biomass
bass_biomass <- subdata[grep("Striped Bass", subdata$variable), ]
bass_biomass <- aggregate(value~Year+Scenario, bass_biomass, sum)

year_id <- seq(1, nrow(amo_unsmooth_lag1), by = 12)

amo0 <- amo_unsmooth_lag0[year_id, ]
amo1 <- amo_unsmooth_lag1[year_id, ]
pcp <- precipitation[year_id, ]
pdsi <- palmer_drought_severity_index[year_id, ]
sst <- kaplan_sst[year_id, ]
effort <- apply(ewe_faa[, grep("menhaden", colnames(ewe_faa))] / ewe_catchability[year_id, grep("menhaden", colnames(ewe_catchability))], 1, sum)
fishing_mortality <- apply(ewe_faa[, grep("menhaden", colnames(ewe_faa))], 1, max)


indices_name <- c("AMO0", "AMO1", "PCP", "PDSI", "SST", "FISHING_EFFORT", "FISHING_MORTALITY", "BASS_BIOMASS")

lm_results <- vector(mode = "list", length = length(subscenarios))

for (i in 1:length(subscenarios)){

  biomass_temp <- subdata[grep("Atlantic Menhaden", subdata$variable), ]
  biomass_temp <- biomass_temp[which(biomass_temp$Scenario == scenario_name[scenario_code %in% subscenarios[i]]), ]

  total_biomass <- aggregate(value~Year, biomass_temp, sum)$value

  temp_data <- data.frame(
    Biomass = total_biomass,
    AMO0 = amo0$scaled_value,
    AMO1 = amo1$scaled_value,
    PCP = pcp$scaled_value, 
    PDSI = pdsi$scaled_value,
    SST = sst$scaled_value,
    FISHING_EFFORT = effort,
    FISHING_MORTALITY = fishing_mortality,
    BASS_BIOMASS = bass_biomass$value[bass_biomass$Scenario == scenario_name[scenario_code %in% subscenarios[i]]]
  )
  lag <- 0:7
  lm_results[[i]] <- matrix(NA, nrow = length(lag), ncol= length(indices_name))
  dimnames(lm_results[[i]]) <- list(paste0("lag", lag), indices_name)
  for (lead_time in 1:length(lag)){
    for (j in 1:length(indices_name)){
    lm_results[[i]][lead_time, j] <- paste0(round(summary(lm(temp_data$Biomass[(lag[lead_time]+1):nrow(temp_data)] ~ temp_data[1:(nrow(temp_data)-lag[lead_time]), j+1]))$coefficients[2, 1], digits = 2), 
                             if(summary(lm(temp_data$Biomass[(lag[lead_time]+1):nrow(temp_data)] ~ temp_data[1:(nrow(temp_data)-lag[lead_time]), j+1]))$coefficients[2, 4] <= 0.05) {"*"})
    }
  }


}

Slope of linear regression model

Biomass of menhaden-like species shows significant positive relationship with AMO1, PDSI, SST, FISHING_EFFORT, FISHING_MORTALITY, and BASS_BIOMASS when lag is 1.

case_id <- 3
knitr::kable(lm_results[[case_id]], caption = paste("Slope of linear regression model from case", scenario_name[case_id]))

{-}

Single-species stock assessment input data simulation {.tabset}

library(RColorBrewer)
library(ggplot2)
library(tidyverse)
library(reshape2)
library(here)
set.seed(999)

options(digits = 2)

# specify working directories ---------------------------------------------

project_path <- here::here()
ewe_output_path <- file.path(project_path, "data", "ewe", "7ages_newsim_final", "ewe7ages_ecosim_final", "ecosim_forcing_pdsi_egg_amo1")

menhadenSA_output_path <- file.path(project_path, "data", "AtlanticMenhadenSA")

# read Atlantic menhaden Beaufort Assessment Model output data --------

menhadenSA_output <- dget(file.path(menhadenSA_output_path, "am019.rdat"))

years <- 1985:2017
ages <- 0:6

# fishery selectivity
fishery_sel <- IFA4EBFM::logistic(
  pattern = "double_logistic",
  x = ages,
  slope_asc = 3.1,
  location_asc = 1.8,
  slope_desc = 0.88,
  location_desc = 0.01
) # CRS 1972 selectivity

# F-at-age
phase1_nyear <- 15
phase2_nyear <- 5
set.seed(9999)
f_full <- c(
  rep(0.15, length=phase1_nyear),
  seq(0.25, 0.7, length = phase2_nyear),
  seq(0.7, 0.2, length = length(years) - phase1_nyear - phase2_nyear)
) *
  exp(rnorm(length(years), mean = 0, sd = 0.2))
fatage <- matrix(NA, nrow = length(f_full), ncol = length(fishery_sel))
rownames(fatage) <- years
colnames(fatage) <- paste0("menhaden", ages)
for (i in 1:length(f_full)) {
  fatage[i, ] <- f_full[i] * fishery_sel
}

write.csv(fatage, file.path(project_path, "data", "ewe", "7ages_newsim_final", "fatage.csv"))

Biomass-at-age

skip_nrows <- 8
species <- 4:10
species_labels <- paste0("age", ages)
species_vec <- paste0("X", species)

# Load biomass data
temp <- scan(file.path(ewe_output_path, "biomass_annual.csv"), what = "", sep = "\n")
data <- temp[-c(1:skip_nrows)]
data <- read.table(
  text = as.character(data),
  sep = ",",
  col.names = read.table(text = temp[skip_nrows], sep = ",")
)

data <- data[, c("year.group", species_vec)]
data[, species_vec] <- data[, species_vec] * 1000000  # biomass in mt
colnames(data) <- c("year", species_labels)

#data[, species_labels] <- data[, species_labels] / rowSums(data[, species_labels])
baa <- reshape2::melt(data, id.vars = c("year"))
colnames(baa) <- c("Year", "Age", "Biomass")

p <- ggplot(data=baa, aes(x=Age, fill=Age, y=Biomass)) +
      geom_bar(position='dodge', stat='identity') +
      facet_wrap(~Year) +
      labs(
        title = "Biomass-at-age",
        x = "Age",
        y = "Proportion"
      )+
      theme_bw()+
      theme(axis.text.x = element_text(angle = 45))
print(p)

Catch

# Create fishery ----------------------------------------------------------

fishery_sample_num <- cbind(
  menhadenSA_output$t.series$acomp.cRs.n[which(menhadenSA_output$t.series$year %in% years)],
  menhadenSA_output$t.series$acomp.cBs.n[which(menhadenSA_output$t.series$year %in% years)],
  menhadenSA_output$t.$acomp.cBn.n[which(menhadenSA_output$t.series$year %in% years)],
  menhadenSA_output$t.series$acomp.cRn.n[which(menhadenSA_output$t.series$year %in% years)]
)
fishery_sample_num[fishery_sample_num==-99999] <- 0

fishery <- create_fishery(
  file_path = file.path(ewe_output_path, "catch_annual.csv"),
  skip_nrows = 8,
  species = 4:10,
  species_labels = paste0("age", ages),
  ewe_years = 0:32,
  data_years = years,
  fleet_num = 1,
  selectivity = NULL,
  CV = rep(0.05, length(years)),
  sample_num = apply(fishery_sample_num, 1, sum),
  waa_path = file.path(ewe_output_path, "weight_annual.csv")
)

# Total catch in biomass
om <- data.frame (
  Year = years, 
  Catch = fishery$om_total_catch_biomass,
  Model = "OM"
)

obs <- data.frame (
  Year = years, 
  Catch = fishery$obs_total_catch_biomass$fleet1,
  Model = "OBS"
)

year_id <- which(menhadenSA_output$t.series$year %in% years)
bam <- data.frame (
  Year = years, 
  Catch = (menhadenSA_output$t.series$L.cBs.ob[year_id] + menhadenSA_output$t.series$L.cBn.ob[year_id] + menhadenSA_output$t.series$L.cRs.ob[year_id] + menhadenSA_output$t.series$L.cRn.ob[year_id]) * 1000,
  Model = "BAM"
)

#data <- rbind(om, obs, bam)
data <- rbind(om, obs)
p <- ggplot(data, aes(x = Year, y = Catch, color = Model)) +
  geom_line() +
  labs(
    title = "Comparison between models",
    x = "Year",
    y = "Total Catch (mt)"
  ) +
  theme_bw()
p

Catch-at-age

om <- fishery$om_caa_abundance/rowSums(fishery$om_caa_abundance)
om$Year <- years
om$Model <- "OM"
om <- reshape2::melt(om, id.vars = c("Year", "Model"))
colnames(om) <- c("Year", "Model", "Age", "Catch")

obs <- fishery$obs_caa_prop$fleet1
obs$Year <- years
obs$Model <- "OBS"
obs <- reshape2::melt(obs, id.vars = c("Year", "Model"))
colnames(obs) <- c("Year", "Model", "Age", "Catch")

data <- list(om=om, obs=obs)
## OM catch-at-age
p <- ggplot(data=data$om, aes(x=Age, fill=Age, y=Catch)) +
      geom_bar(position='dodge', stat='identity') +
      facet_wrap(~Year) +
      labs(
        title = "OM catch-at-age",
        x = "Age",
        y = "Proportion"
      )+
      theme_bw()+
      theme(axis.text.x = element_text(angle = 45)) 

p 

## Observed catch-at-age
p <- ggplot(data=data$obs, aes(x=Age, fill=Age, y=Catch)) +
      geom_bar(position='dodge', stat='identity') +
      facet_wrap(~Year) +
      labs(
        title = "Observed catch-at-age",
        x = "Age",
        y = "Proportion"
      )+
      theme_bw()+
      theme(axis.text.x = element_text(angle = 45))

p

# Create survey ----------------------------------------------------------

# selectivity settings
survey_num <- 2

survey_name <- c("survey1", "survey2")

# set up survey time
# Need to check Table 26 from BAM assessment: Length cutoffs used to distinguish age-0 from age-1+ Atlantic menhaden at different regions.

survey1_year <- 1990:2017 # Adult Index (survey1): age 1+ fish; September - January; Time of the year when menhaden were most abundant in this region

survey2_year <- 1985:2017 # Adult Index (survey2): age 1+ fish; March - May

survey_time <- list(
  survey1 = data.frame(
    year = survey1_year,
    month = rep(10, length(survey1_year)) # Oct 15
  ),
  survey2 = data.frame(
    year = survey2_year,
    month = rep(4, length(survey2_year)) # April 15
  )
)

# set up survey selectivity
# survey1_sel <- IFA4EBFM::logistic(
#   pattern = "double_logistic",
#   x = ages,
#   slope_asc = 2.2,
#   location_asc = 3.0,
#   slope_desc = 3.0,
#   location_desc = 3.5
# )
# 
# survey1_sel <- IFA4EBFM::logistic(
#   pattern = "simple_logistic",
#   x = ages,
#   slope_asc = 2.2,
#   location_asc = 3.0
# )

survey1_sel <- IFA4EBFM::logistic(
  pattern = "double_logistic",
  x = ages,
  slope_asc = 4.3,
  location_asc = 2.3,
  slope_desc = 3.5,
  location_desc = 2.3
)

survey2_sel <- IFA4EBFM::logistic(
  pattern = "double_logistic",
  x = ages,
  slope_asc = 4.3,
  location_asc = 2.3,
  slope_desc = 3.5,
  location_desc = 2.3
)


survey_selectivity <- list(
  survey1 = as.data.frame(
    matrix(rep(survey1_sel, times = length(years)), ncol = length(ages), byrow = TRUE),
    row.names = years
  ),
  survey2 = as.data.frame(
    matrix(rep(survey2_sel, times = length(years)), ncol = length(ages), byrow = TRUE),
    row.names = years,
  )
)

survey_selectivity <- lapply(survey_selectivity, setNames, paste("age", ages))

# set up catchability
survey_catchability <- list(
  survey1 = menhadenSA_output$t.series$q.nad[which(menhadenSA_output$t.series$year %in% years)],
  survey2 = menhadenSA_output$t.series$q.mad[which(menhadenSA_output$t.series$year %in% years)]
)

survey_catchability <- lapply(survey_catchability, setNames, years)

# set up CV
# survey_CV <- list(
#   survey1 = menhadenSA_output$t.series$cv.U.nad[which(menhadenSA_output$t.series$year %in% years)],
#   survey2 = menhadenSA_output$t.series$cv.U.mad[which(menhadenSA_output$t.series$year %in% years)]
# )

survey_CV <- list(
  survey1 = rep(0.1, length(years)),
  survey2 = rep(0.1, length(years))
)

survey_CV <- lapply(survey_CV, setNames, years)

# set up sample number
survey_sample_num <- list(
  survey1 = menhadenSA_output$t.series$lcomp.nad.nfish[which(menhadenSA_output$t.series$year %in% years)],
  #survey2 = menhadenSA_output$t.series$lcomp.mad.nfish[which(menhadenSA_output$t.series$year %in% years)]
  survey2 = rep(1000, length= length(years))
)
survey_sample_num <- lapply(survey_sample_num, setNames, years)

for (i in 1:length(survey_sample_num)) {
  survey_sample_num[[i]][survey_sample_num[[i]] == -99999] <- NA
}

# set up age-length population structure
length_bin <- seq(10.0, 600, 10) / 10 # in cm
mid_length_bin <- seq(15, 605, 10) / 10 # in cm
nbin <- length(length_bin)
bin_width <- 1

length_CV <- list(
  survey1 = 0.12,
  survey2 = 0.17
)

# Create survey
survey <- IFA4EBFM::create_survey(
  file_path = file.path(ewe_output_path, "biomass_monthly.csv"),
  skip_nrows = 8,
  species = 4:10,
  species_labels = paste0("age", ages),
  years = years,
  survey_num = survey_num,
  survey_time = survey_time,
  selectivity = survey_selectivity,
  catchability = survey_catchability,
  CV = survey_CV,
  sample_num = survey_sample_num,
  waa_path = file.path(ewe_output_path, "weight_monthly.csv"),
  length_bin = length_bin,
  mid_length_bin = mid_length_bin,
  nbin = nbin,
  bin_width = bin_width,
  length_CV = length_CV
)

Survey selectivity patterns

plot(0:6, as.numeric(survey_selectivity$survey1[1, ]), 
     type="o", pch=1, col=1, lty=1,
     xlab = "Age", ylab = "Selectivity")
lines(0:6, as.numeric(survey_selectivity$survey2[1, ]), pch=2, col=2, lty=2, type="o")
legend("topright",
       names(survey_selectivity),
       pch= 1:2,
       col = 1:2, 
       lty=1:2,
       bty="n"
       )

Survey index over time

om <- data.frame(matrix(ncol=4, nrow=0))
for (i in 1:length(survey$om_abundance_index)){
  temp <- data.frame(
  Year = as.numeric(names(survey$om_abundance_index[[i]])), 
  Index = survey$om_abundance_index[[i]],
  Model = "OM",
  Survey = names(survey$om_abundance_index)[i]
)
  om <- rbind(om, temp)
}

obs <- data.frame(matrix(ncol=4, nrow=0))
for (i in 1:length(survey$obs_abundance_index)){
  temp <- data.frame(
  Year = as.numeric(names(survey$obs_abundance_index[[i]])), 
  Index = survey$obs_abundance_index[[i]],
  Model = "OBS",
  Survey = names(survey$obs_abundance_index)[i]
)
  obs <- rbind(obs, temp)
}

data <- rbind(om, obs)
p <- ggplot(data, aes(x = Year, y = Index, color = Model)) +
  geom_line() +
  facet_wrap(~Survey, scales = "free_y") +
  labs(
    title = "Comparison between models",
    x = "Year",
    y = "Abundance Index"
  ) +
  theme_bw()
p

Survey length composition

## Survey 1
obs <- as.data.frame(survey$obs_lencomp_proportion_ss3[[1]])
obs$Year <- as.numeric(row.names(survey$obs_lencomp_proportion_ss3[[1]]))
obs$Model <- "OBS"
obs$Survey <- names(survey$obs_lencomp_proportion_ss3)[1]

obs <- na.omit(obs)

obs <- reshape2::melt(obs, id.vars = c("Year", "Model", "Survey"))
colnames(obs) <- c("Year", "Model", "Survey", "Length_Bin", "Value")

obs$Length_Bin <- as.numeric(as.character(obs$Length_Bin))
p <- ggplot(data=obs, aes(x=Length_Bin, fill=Length_Bin, y=Value)) +
  geom_bar(position='dodge', stat='identity') +
  facet_wrap(~Year) +
  labs(
    title = names(survey$obs_lencomp_proportion_ss3)[1],
    x = "Length Bins (cm)",
    y = "Proportion"
  )+
  theme_bw()
p

## Survey 2
obs <- as.data.frame(survey$obs_lencomp_proportion_ss3[[2]])
obs$Year <- as.numeric(row.names(survey$obs_lencomp_proportion_ss3[[2]]))
obs$Model <- "OBS"
obs$Survey <- names(survey$obs_lencomp_proportion_ss3)[2]

obs <- na.omit(obs)

obs <- reshape2::melt(obs, id.vars = c("Year", "Model", "Survey"))
colnames(obs) <- c("Year", "Model", "Survey", "Length_Bin", "Value")

obs$Length_Bin <- as.numeric(as.character(obs$Length_Bin))
p <- ggplot(data=obs, aes(x=Length_Bin, fill=Length_Bin, y=Value)) +
  geom_bar(position='dodge', stat='identity') +
  facet_wrap(~Year) +
  labs(
    title = names(survey$obs_lencomp_proportion_ss3)[2],
    x = "Length Bins (cm)",
    y = "Proportion"
  )+
  theme_bw()
p

biodata <- create_biodata(nsex=1, narea=1, ages=ages, years=years,
                          length_bin=length_bin, mid_length_bin=mid_length_bin,
                          nbin=nbin, bin_width=bin_width, length_CV=length_CV,
                          annual_weight_path=file.path(ewe_output_path, "weight_annual.csv"),
                          monthly_weight_path=file.path(ewe_output_path, "weight_monthly.csv"),
                          species=4:10,
                          species_labels=paste0("age", ages),
                          skip_nrows=8,
                          lw_a=0.01, lw_b=3,
                          k=0.331,
                          t0 = -0.1,
                          winf = 0.237,
                          maturity_at_age=c(0.0, 0.1, 0.5, 0.9, 1.0, 1.0, 1.0), # From both BAM and EwE
                          natural_mortality_at_age=c(1.76, 1.31, 1.03, 0.9, 0.81, 0.76, 0.72) # From both BAM and EwE
)

sa_data <- list(
    fishery = fishery,
    survey = survey,
    biodata = biodata
)

{-}

Stock Synthesis 3 results {.tabset}

# Install required R packages -------------------------------------
devtools::load_all()
required_pkg <- c(
  "devtools", "here"
)

pkg_to_install <- required_pkg[!(required_pkg %in%
  installed.packages()[, "Package"])]
if (length(pkg_to_install)) install.packages(pkg_to_install)

lapply(required_pkg, library, character.only = TRUE)

devtools::install_github("r4ss/r4ss",
  ref = "c53f82fcfb3f54296d79ba3a4163990150981285"
)
library(r4ss)

# Case 0: default stock assessment run ----------------------------
# Load simulated input data
file_path <- here::here("data", "ewe", "7ages_newsim_final", "ewe7ages_ecosim_final", "ecosim_forcing_pdsi_egg_amo1")
source(here::here("Rscript", "load_indices.R"))
file_path <- here::here("data", "data_rich", "OM2")
originalwd <- getwd()

depletion <- FALSE
depletion_ratio <- 0.8 # tried values between 0.3 - 0.8

model_year <- 1985:2012
projection_year <- 2013:2017

generate_ss3(
  file_path = file_path,
  r0 = 12, steepness = 0.5, sigmar = 0.25,
  projection_f = 0.75, projection_catch = NULL,
  sa_data = sa_data, model_year = model_year,
  projection_year = projection_year,
  use_depletion = FALSE, depletion_ratio = NULL, 
  initial_equilibrium_catch = TRUE
)

# Run SS3 ----------------------------------------------------------

setwd(file_path)
system(paste(file.path(file_path, "ss_win.exe"), file.path(file_path, "data.ss"), sep = " "), show.output.on.console = TRUE)

# get FMSY
ss3list <- SS_output(
  dir = file_path,
  verbose = TRUE,
  printstats = TRUE
)

fmsy <- ss3list$derived_quants$Value[ss3list$derived_quants$Label == "annF_MSY"]

# Re-run SS3 with projection_f = fmsy
generate_ss3(
  file_path = file_path,
  r0 = 12, steepness = 0.5, sigmar = 0.25,
  projection_f = fmsy, projection_catch = NULL,
  sa_data = sa_data, model_year = model_year,
  projection_year = projection_year,
  use_depletion = FALSE, depletion_ratio = NULL,
  initial_equilibrium_catch = TRUE
)

system(paste(file.path(file_path, "ss_win.exe"), file.path(file_path, "data.ss"), sep = " "), show.output.on.console = TRUE)
on.exit(setwd(originalwd), add = TRUE)
# plot r4ss -------------------------------------------------------

# read the model output and print diagnostic messages
ss3list <- SS_output(
  dir = file_path,
  verbose = TRUE,
  printstats = TRUE
)

# plots the results
# SS_plots(ss3list)

Compare "true" values from OM and SS3 results

# Compare SS3 estimate with EwE "true" values ---------------------
# Comparisons
functional_groups <- c(
  "StripedBass0",
  "StripedBass2_5",
  "StripedBass6",
  "AtlanticMenhaden0",
  "AtlanticMenhaden1",
  "AtlanticMenhaden2",
  "AtlanticMenhaden3",
  "AtlanticMenhaden4",
  "AtlanticMenhaden5",
  "AtlanticMenhaden6",
  "SpinyDogfish",
  "BluefishJuvenile",
  "BluefishAdult",
  "WeakfishJuvenile",
  "WeakfishAdult",
  "AtlanticHerring0_1",
  "AtlanticHerring2",
  "Anchovies",
  "Benthos",
  "Zooplankton",
  "Phytoplankton",
  "Detritus"
)

age_name <- paste0("AtlanticMenhaden", 0:6)

# Biomass
biomass <- read_ewe_output(
  file_path = ewe_output_path,
  file_names = "biomass_monthly.csv",
  skip_nrows = 8,
  plot = FALSE,
  figure_titles = NULL,
  functional_groups = functional_groups,
  figure_colors = NULL
)

time_id <- seq(1, nrow(biomass[[1]]), by = 12)
biomass_ewe <- apply(biomass[[1]][, age_name], 1, sum) * 1000000

par(mfrow = c(2, 2))
model_year_id <- which(ss3list$timeseries$Yr %in% model_year)
projection_year_id <- which(ss3list$timeseries$Yr %in% projection_year)
ylim <- range(biomass_ewe[time_id], ss3list$timeseries$Bio_all)
plot(c(model_year, projection_year),
  biomass_ewe[time_id],
  xlab = "Year", ylab = "Biomass (mt)",
  ylim = ylim,
  pch = 16
)
lines(model_year, ss3list$timeseries$Bio_all[model_year_id], lty = 1)
#par(new = TRUE)
#plot(model_year, ss3list$timeseries$Bio_all[model_year_id], lty = 1, type = "l", axes = FALSE, xlab = "", ylab = "")
lines(projection_year,
  ss3list$timeseries$Bio_all[projection_year_id],
  lty = 2
)
#axis(side = 4, at = pretty(range(ss3list$timeseries$Bio_all[model_year_id])))
legend("topleft",
  c("EWE", "SS3 Estimation", "SS3 projection"),
  bty = "n",
  pch = c(16, NA, NA),
  lty = c(NA, 1, 2)
)
cor(biomass_ewe[time_id], c(ss3list$timeseries$Bio_all[model_year_id], ss3list$timeseries$Bio_all[projection_year_id]))

# Abundance
waa <- read_ewe_output(
  file_path = ewe_output_path,
  file_names = "weight_monthly.csv",
  skip_nrows = 8,
  plot = FALSE,
  figure_titles = NULL,
  functional_groups = functional_groups,
  figure_colors = NULL
)

abundance_ewe <- apply(biomass[[1]][, age_name] * 1000000 / waa[[1]][, age_name], 1, sum)[time_id]

ylim <- range(abundance_ewe, ss3list$timeseries$`SmryNum_SX:1_GP:1` * 1000)
plot(c(model_year, projection_year),
  abundance_ewe,
  xlab = "Year", ylab = "Abundance",
  ylim = ylim,
  pch = 16
)
lines(model_year, ss3list$timeseries$`SmryNum_SX:1_GP:1`[model_year_id] * 1000, lty = 1)
#par(new = TRUE)
#plot(model_year, ss3list$timeseries$`SmryNum_SX:1_GP:1`[model_year_id] * 1000, lty = 1, type = "l", axes = FALSE, xlab = "", ylab = "")
#axis(side = 4, at = pretty(range(ss3list$timeseries$`SmryNum_SX:1_GP:1`[model_year_id] * 1000)))
lines(projection_year,
  ss3list$timeseries$`SmryNum_SX:1_GP:1`[projection_year_id] * 1000,
  lty = 2
)
legend("topright",
  c("EWE", "SS3 Estimation", "SS3 projection"),
  bty = "n",
  pch = c(16, NA, NA),
  lty = c(NA, 1, 2)
)

cor(abundance_ewe, c(ss3list$timeseries$`SmryNum_SX:1_GP:1`[model_year_id] * 1000, ss3list$timeseries$`SmryNum_SX:1_GP:1`[projection_year_id] * 1000))
# Recruitment
waa <- read_ewe_output(
  file_path = ewe_output_path,
  file_names = "weight_monthly.csv",
  skip_nrows = 8,
  plot = FALSE,
  figure_titles = NULL,
  functional_groups = functional_groups,
  figure_colors = NULL
)

recruit_ewe <- biomass[[1]][, "AtlanticMenhaden0"][time_id] * 1000000 / waa[[1]][, "AtlanticMenhaden0"][time_id]

recruit_ss3 <- data.frame(
  year = c(model_year, projection_year),
  recruit = ss3list$natage$`1`[which(ss3list$natage$`Beg/Mid` == "B" & ss3list$natage$Yr %in% c(model_year, projection_year))] * 1000
)


ylim <- range(recruit_ewe, recruit_ss3$recruit)
plot(c(model_year, projection_year),
  recruit_ewe,
  xlab = "Year", ylab = "Abundance of age-0 fish",
  ylim = ylim,
  pch = 16
)
lines(model_year, recruit_ss3$recruit[recruit_ss3$year %in% model_year],
     lty = 1)
# par(new = TRUE)
# plot(model_year, recruit_ss3$recruit[recruit_ss3$year %in% model_year],
#   lty = 1, type = "l", axes = FALSE, xlab = "", ylab = ""
# )
# axis(side = 4, at = pretty(range(recruit_ss3$recruit[recruit_ss3$year %in% model_year])))
lines(projection_year,
  recruit_ss3$recruit[recruit_ss3$year %in% projection_year],
  lty = 2
)
legend("topright",
  c("EWE", "SS3 Estimation", "SS3 projection"),
  bty = "n",
  pch = c(16, NA, NA),
  lty = c(NA, 1, 2)
)

# Mortality
# mortality <- read.csv(file = here::here("data", "ewe", "7ages", "updated_f.csv"))
# mortality_ewe <- apply(mortality[seq(1, nrow(mortality), by=12), grep("AtlanticMenhaden", colnames(mortality))], 1, max)

mortality <- read.csv(file = here::here("data", "ewe", "7ages_newsim_final", "fatage.csv"))
mortality_ewe <- apply(mortality[, 2:ncol(mortality)], 1, max)
# mortality_ewe <- apply(mortality[, 2:ncol(mortality)], 1, mean)

ylim <- range(mortality_ewe, ss3list$timeseries$`F:_1`)
plot(c(model_year, projection_year),
  mortality_ewe,
  xlab = "Year", ylab = "Mortality",
  ylim = ylim,
  pch = 16
)
lines(model_year, ss3list$timeseries$`F:_1`[model_year_id],
     lty = 1)
# par(new = TRUE)
# plot(model_year, ss3list$timeseries$`F:_1`[model_year_id],
#   lty = 1, type = "l", , axes = FALSE, xlab = "", ylab = ""
# )
# axis(side = 4, at = pretty(range(ss3list$timeseries$`F:_1`[model_year_id])))
lines(projection_year,
  ss3list$timeseries$`F:_1`[projection_year_id],
  lty = 2
)
legend("topright",
  c("EWE F", "SS3 F", "SS3 FMSY"),
  bty = "n",
  pch = c(16, NA, NA),
  lty = c(NA, 1, 2)
)

Linear regression models

# Linear regression models from cases 1 - 4 using "true" values from EwE------------------------

# Biomass
biomass <- read_ewe_output(
  file_path = ewe_output_path,
  file_names = "biomass_monthly.csv",
  skip_nrows = 8,
  plot = FALSE,
  figure_titles = NULL,
  functional_groups = functional_groups,
  figure_colors = NULL
)

time_id <- seq(1, nrow(biomass[[1]]), by = 12)
biomass_ewe <- apply(biomass[[1]][, age_name], 1, sum) 


lm_slope <- data.frame(
  case = "True indices",
  amo = NA,
  pdsi = NA,
  bassB = NA,
  effort =NA
)


year_id <- seq(1, nrow(amo_unsmooth_lag1), by = 12)
index_year <- c(model_year, projection_year)
lag = 1
biomass_lag_id <- (1+lag):length(index_year)
index_lag_id <- 1:(length(index_year)-lag)

amo <- amo_unsmooth_lag1[year_id, ]
amo_lm <- lm(biomass_ewe[time_id][biomass_lag_id] ~ amo$scaled_value[index_lag_id])
amo_fit <- fitted(amo_lm)
lm_slope$amo <- paste0(
  round(summary(amo_lm)$coefficients[2, 1], digits = 2),
  if (summary(amo_lm)$coefficients[2, 4] <= 0.05) {
    "*"
  }
)

pdsi <- palmer_drought_severity_index[year_id, ]
pdsi_lm <- lm(biomass_ewe[time_id][biomass_lag_id] ~ pdsi$scaled_value[index_lag_id])
pdsi_fit <- fitted(pdsi_lm)
lm_slope$pdsi <- paste0(
  round(summary(pdsi_lm)$coefficients[2, 1], digits = 2),
  if (summary(pdsi_lm)$coefficients[2, 4] <= 0.05) {
    "*"
  }
)

bassB <- bass_bio[bass_bio$Year %in% year_id, ]
bassB_lm <- lm(biomass_ewe[time_id][biomass_lag_id] ~ bassB$bass_bio[index_lag_id])
bassB_fit <- fitted(bassB_lm)
lm_slope$bassB <- paste0(
  round(summary(bassB_lm)$coefficients[2, 1], digits = 2),
  if (summary(bassB_lm)$coefficients[2, 4] <= 0.05) {
    "*"
  }
)

effort <- apply(ewe_faa[, grep("menhaden", colnames(ewe_faa))] / ewe_catchability[year_id, grep("menhaden", colnames(ewe_catchability))], 1, sum)
effort_lm <- lm(biomass_ewe[time_id][biomass_lag_id] ~ effort[index_lag_id])
effort_fit <- fitted(effort_lm)
lm_slope$effort <- paste0(
  round(summary(effort_lm)$coefficients[2, 1], digits = 2),
  if (summary(effort_lm)$coefficients[2, 4] <= 0.05) {
    "*"
  }
)


par(mfrow = c(2, 2))

plot(amo$scaled_value[index_lag_id], biomass_ewe[time_id][biomass_lag_id],
  xlab = "AMO values",
  ylab = "Menhaden-like species B from OM"
)
lines(amo$scaled_value[index_lag_id], amo_fit, lty = 2, col = "blue")

plot(pdsi$scaled_value[index_lag_id], biomass_ewe[time_id][biomass_lag_id],
  xlab = "PDSI values",
  ylab = "Menhaden-like species B from OM"
)
lines(pdsi$scaled_value[index_lag_id], pdsi_fit, lty = 2, col = "blue")

plot(bassB$bass_bio[index_lag_id], biomass_ewe[time_id][biomass_lag_id],
  xlab = "Biomass of Striped bass",
  ylab = "Menhaden-like species B from OM"
)
lines(bassB$bass_bio[index_lag_id], bassB_fit, lty = 2, col = "blue")

plot(effort[index_lag_id], biomass_ewe[time_id][biomass_lag_id],
  xlab = "Menhaden fishing effort",
  ylab = "Menhaden-like species B from OM"
)
lines(effort[index_lag_id], effort_fit, lty = 2, col = "blue")

knitr::kable(lm_slope)

# Projection: cases 1 - 4 ---------------------------

lm_slope <- data.frame(
  case = "PDSI+AMO1 Case",
  projection_year = 1:length(projection_year),
  amo = NA,
  pdsi = NA,
  bassB = NA,
  effort = NA
)

for (projection_year_id in 1:length(projection_year)) {

  # abundance
  # if (projection_year_id == 1) {
  #   menhaden_b <- ss3list$timeseries$`SmryNum_SX:1_GP:1`[model_year_id]
  # } else {
  #   menhaden_b <- c(ss3list$timeseries$`SmryNum_SX:1_GP:1`[model_year_id], ss3list$timeseries$`SmryNum_SX:1_GP:1`[1:(projection_year_id - 1)])
  # }

  # biomass
  if (projection_year_id ==1) {
    menhaden_b <- ss3list$timeseries$Bio_all[model_year_id]
  } else {
    menhaden_b <- c(ss3list$timeseries$Bio_all[model_year_id], ss3list$timeseries$Bio_all[1:(projection_year_id-1)])
  }


  year_id <- seq(1, nrow(amo_unsmooth_lag1), by = 12)[1:(length(model_year) + projection_year_id - 1)]

  if (projection_year_id == 1) {
    index_year <- model_year
  } else {
    index_year <- c(model_year, projection_year[1:(projection_year_id - 1)])
  }


  # Linear regression model -----------------------------------------

  lag = 1
  biomass_lag_id <- (1+lag):length(index_year)
  index_lag_id <- 1:(length(index_year)-lag)

  amo <- amo_unsmooth_lag1[year_id, ]
  amo_lm <- lm(menhaden_b[biomass_lag_id] ~ amo$scaled_value[index_lag_id])
  amo_fit <- fitted(amo_lm)
  lm_slope$amo <- paste0(
    round(summary(amo_lm)$coefficients[2, 1], digits = 2),
    if (summary(amo_lm)$coefficients[2, 4] <= 0.05) {
      "*"
    }
  )

  pdsi <- palmer_drought_severity_index[year_id, ]
  pdsi_lm <- lm(menhaden_b[biomass_lag_id] ~ pdsi$scaled_value[index_lag_id])
  pdsi_fit <- fitted(pdsi_lm)
  lm_slope$pdsi <- paste0(
    round(summary(pdsi_lm)$coefficients[2, 1], digits = 2),
    if (summary(pdsi_lm)$coefficients[2, 4] <= 0.05) {
      "*"
    }
  )

  bassB <- bass_bio[bass_bio$Year %in% year_id, ]
  bassB_lm <- lm(menhaden_b[biomass_lag_id] ~ bassB$bass_bio[index_lag_id])
  bassB_fit <- fitted(bassB_lm)
  lm_slope$bassB <- paste0(
    round(summary(bassB_lm)$coefficients[2, 1], digits = 2),
    if (summary(bassB_lm)$coefficients[2, 4] <= 0.05) {
      "*"
    }
  )

  effort <- apply(ewe_faa[ewe_faa$X %in% ewe_faa$X[1:length(year_id)] , grep("menhaden", colnames(ewe_faa))] / ewe_catchability[year_id, grep("menhaden", colnames(ewe_catchability))], 1, sum)
  effort_lm <- lm(menhaden_b[biomass_lag_id] ~ effort[index_lag_id])
  effort_fit <- fitted(effort_lm)
  lm_slope$effort <- paste0(
    round(summary(effort_lm)$coefficients[2, 1], digits = 2),
    if (summary(effort_lm)$coefficients[2, 4] <= 0.05) {
      "*"
    }
  )

  if (projection_year_id == length(projection_year)) {
    par(mfrow = c(2, 2))

    plot(amo$scaled_value[index_lag_id], menhaden_b[biomass_lag_id],
      xlab = "AMO values",
      ylab = "Biomass of menhaden-like species"
    )
    lines(amo$scaled_value[index_lag_id], amo_fit, lty = 2, col = "blue")

    plot(pdsi$scaled_value[index_lag_id], menhaden_b[biomass_lag_id],
      xlab = "PDSI values",
      ylab = "Biomass of menhaden-like species"
    )
    lines(pdsi$scaled_value[index_lag_id], pdsi_fit, lty = 2, col = "blue")

    plot(bassB$bass_bio[index_lag_id], menhaden_b[biomass_lag_id],
      xlab = "Biomass of Striped bass",
      ylab = "Biomass of menhaden-like species"
    )
    lines(bassB$bass_bio[index_lag_id], bassB_fit, lty = 2, col = "blue")

    plot(effort[index_lag_id], menhaden_b[biomass_lag_id],
      xlab = "Menhaden fishing effort",
      ylab = "Biomass of menhaden-like species"
    )
    lines(effort[index_lag_id], effort_fit, lty = 2, col = "blue")
  }


  # status of indicators --------------------------------------------

  amo_soi <- calc_soi(
    indicator_data = amo$scaled_value,
    slope = coef(amo_lm)[2]
  )

  pdsi_soi <- calc_soi(
    indicator_data = pdsi$scaled_value,
    slope = coef(pdsi_lm)[2]
  )

  bassB_soi <- calc_soi(
    indicator_data = bassB$bass_bio,
    slope = coef(bassB_lm)[2]
  )

  effort_soi <- calc_soi(
    indicator_data = effort,
    slope = coef(effort_lm)[2]
  )

  if (projection_year_id == 1) {
    scaled_data <- data.frame(
      year = model_year,
      projection_year_id = projection_year_id,
      amo = scale(amo$scaled_value)[, 1],
      pdsi = scale(pdsi$scaled_value)[, 1],
      bassB = scale(bassB$bass_bio)[, 1],
      effort = scale(effort)[, 1],
      menhadenB = scale(menhaden_b)[, 1]
    )

    soi_data <- data.frame(
      year = model_year,
      projection_year_id = projection_year_id,
      amo = amo_soi,
      pdsi = pdsi_soi,
      bass_b = bassB_soi,
      effort = effort_soi
    )
  } else {
    scaled_data <- rbind(
      scaled_data,
      data.frame(
        year = index_year,
        projection_year_id = projection_year_id,
        amo = scale(amo$scaled_value)[, 1],
        pdsi = scale(pdsi$scaled_value)[, 1],
        bassB = scale(bassB$bass_bio)[, 1],
        effort = scale(effort)[, 1],
        menhadenB = scale(menhaden_b)[, 1]
      )
    )

    soi_data <- rbind(
      soi_data,
      data.frame(
        year = index_year,
        projection_year_id = projection_year_id,
        amo = amo_soi,
        pdsi = pdsi_soi,
        bass_b = bassB_soi,
        effort = effort_soi
      )
    )
  }



  # Adjusted FMSY ----------------------------------------------------

  if (projection_year_id == 1) {
    Bt_BMSY <- ss3list$Kobe$B.Bmsy[length(model_year)]
  } else {
    Bt_BMSY <- ss3list$Kobe$B.Bmsy[length(model_year) + projection_year_id - 1]
  }

  ss3_fmsy <- rnorm(1000, ss3list$derived_quants$Value[ss3list$derived_quants$Label == "annF_MSY"], 0.5)
  ss3_Bt_BMSY <- rep(Bt_BMSY, 1000)

  ss3_fmsy[ss3_fmsy < 0] <- 0.1
  amo_fmsy <- adjust_projection_jabba(
    FMSY = ss3_fmsy,
    soi = tail(amo_soi, n = 1),
    Bt_BMSY = ss3_Bt_BMSY
  )
  pdsi_fmsy <- adjust_projection_jabba(
    FMSY = ss3_fmsy,
    soi = tail(pdsi_soi, n = 1),
    Bt_BMSY = ss3_Bt_BMSY
  )
  bassB_fmsy <- adjust_projection_jabba(
    FMSY = ss3_fmsy,
    soi = tail(bassB_soi, n = 1),
    Bt_BMSY = ss3_Bt_BMSY
  )
  effort_fmsy <- adjust_projection_jabba(
    FMSY = ss3_fmsy,
    soi = tail(effort_soi, n = 1),
    Bt_BMSY = ss3_Bt_BMSY
  )

  if (projection_year_id == 1) {
    fmsy_data <- data.frame(
      iter = 1:length(amo_fmsy),
      projection_year_id = projection_year_id,
      ss3 = ss3_fmsy,
      amo = amo_fmsy,
      pdsi = pdsi_fmsy,
      bassB = bassB_fmsy,
      effort = effort_fmsy
    )
  } else {
    fmsy_data <- rbind(
      fmsy_data,
      data.frame(
        iter = 1:length(amo_fmsy),
        projection_year_id = projection_year_id,
        ss3 = ss3_fmsy,
        amo = amo_fmsy,
        pdsi = pdsi_fmsy,
        bassB = bassB_fmsy,
        effort = effort_fmsy
      )
    )
  }
}



scaled_data_melt <- reshape2::melt(
  scaled_data,
  id = c("year", "projection_year_id")
)
scaled_data_melt$projection_year_id <- as.factor(scaled_data_melt$projection_year_id)

soi_data_melt <- reshape2::melt(
  soi_data,
  id = c("year", "projection_year_id")
)
soi_data_melt$projection_year_id <- as.factor(soi_data_melt$projection_year_id)

F adjustment

fmsy_data_melt <- reshape2::melt(
  fmsy_data,
  id = c("iter", "projection_year_id")
)
fmsy_data_melt$projection_year_id <- as.factor(fmsy_data_melt$projection_year_id)

fmsy_data_melt_median <- aggregate(value ~ projection_year_id + variable, data = fmsy_data_melt, median)
fmsy_data_melt_median$projection_year_id <- as.numeric(fmsy_data_melt_median$projection_year_id)

p <- ggplot(fmsy_data_melt, aes(x = variable, y = value, color = projection_year_id)) +
  geom_boxplot() +
  labs(
    title = "FMSY from SS3 and adjusted FMSY",
    x = "",
    y = "FMSY"
  ) +
  theme_bw()

p

{-}

Bai-Li-NOAA/IFA4EBFM documentation built on May 1, 2024, 9:24 p.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

Bai-Li-NOAA/IFA4EBFM
Interdisciplinary Forecasting Approach for Ecosystem-Based Fisheries Management

In Bai-Li-NOAA/IFA4EBFM: Interdisciplinary Forecasting Approach for Ecosystem-Based Fisheries Management

Simulation scenarios {.tabset}

Environmental drivers

{-}

Biomass over time {.tabset}

Biomass trends from cases that use forcing function and AMO index through egg production

Total biomass of menhaden-like species from cases that use forcing function and AMO index through egg production

Slope of linear regression model

{-}

Single-species stock assessment input data simulation {.tabset}

Biomass-at-age

Catch

Catch-at-age

Survey selectivity patterns

Survey index over time

Survey length composition

{-}

Stock Synthesis 3 results {.tabset}

Compare "true" values from OM and SS3 results

Linear regression models

F adjustment

{-}

R Package Documentation

Browse R Packages

We want your feedback!

Bai-Li-NOAA/IFA4EBFM Interdisciplinary Forecasting Approach for Ecosystem-Based Fisheries Management

In Bai-Li-NOAA/IFA4EBFM: Interdisciplinary Forecasting Approach for Ecosystem-Based Fisheries Management

Simulation scenarios {.tabset}

Environmental drivers

{-}

Biomass over time {.tabset}

Biomass trends from cases that use forcing function and AMO index through egg production

Total biomass of menhaden-like species from cases that use forcing function and AMO index through egg production

Slope of linear regression model

{-}

Single-species stock assessment input data simulation {.tabset}

Biomass-at-age

Catch

Catch-at-age

Survey selectivity patterns

Survey index over time

Survey length composition

{-}

Stock Synthesis 3 results {.tabset}

Compare "true" values from OM and SS3 results

Linear regression models

F adjustment

{-}

R Package Documentation

Browse R Packages

We want your feedback!

Bai-Li-NOAA/IFA4EBFM
Interdisciplinary Forecasting Approach for Ecosystem-Based Fisheries Management