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

Case 0: stock assessment base run (terminal year = 2008)

# Load simulated input data ---------------------------------------------
source(here::here("Rscript", "simulationOM3.R")) # OM reference points
maturity_matrix <- sa_data$biodata$maturity_matrix

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"
)
ages <- 0:6
age_name <- paste0("AtlanticMenhaden", ages)
start_year <- 1985
projection_nyear <- 5
model_year <- start_year:terminal_year
projection_year <- (terminal_year + 1):(terminal_year + projection_nyear)

figure_path <- here::here("figure", "data_moderate")
if (!dir.exists(figure_path)) dir.create(figure_path)

# Reference points scenarios
# Maximum relative F = 7
# No. of steps = 700
# No. of trial years = 100
# Option 1: Search by fleets or groups
# Option 2: Full compensation or stationary system
if (add_environmental_effects == TRUE){
  msy_file_path <- here::here("data", "ewe", "7ages_ecosim_om_msy", "msy_forcing_pdsi_egg_amo1_fleet_details")
} else {
  msy_file_path <- here::here("data", "ewe", "7ages_ecosim_om_msy", "msy_base_run_fleet_details")
}

compensation_msy <- read_ewe_reference_points_vec(
  file_path = msy_file_path,
  file_names = c(
    "F.menhaden_Catch_FullCompensation.csv",
    "F.menhaden_B_FullCompensation.csv"
  ),
  reference_points_scenario = "compensation",
  functional_groups = functional_groups,
  key_functional_group = "AtlanticMenhaden",
  ages = ages,
  maturity_vec = maturity_matrix[1,]
)

stationary_msy <- read_ewe_reference_points_vec(
  file_path = msy_file_path,
  file_names = c(
    "F.menhaden_Catch_StationarySystem.csv",
    "F.menhaden_B_StationarySystem.csv"
  ),
  reference_points_scenario = "compensation",
  functional_groups = functional_groups,
  key_functional_group = "AtlanticMenhaden",
  ages = ages,
  maturity_vec = maturity_matrix[1,]
)

if (add_environmental_effects == TRUE){
  fbase_file_path <- here::here("data", "ewe", "7ages_ecosim_om_msy", "msy_forcing_pdsi_egg_amo1_fleet")
} else {
  fbase_file_path <- here::here("data", "ewe", "7ages_ecosim_om_msy", "msy_base_run_fleet")
}
temp <- scan(file.path(fbase_file_path, "FMSY_FullCompensation.csv"), what = "", sep = "\n")
fmsy_data <- temp[-c(1:9)]
fmsy_data <- read.table(
    text = as.character(fmsy_data),
    sep = ",", 
    col.names = c(
      "Group", "TL", "Fbase", "Cbase", "Vbase", 
      "FmsyFound", "Fmsy", "Cmsy", "Vmsy", "CmsyAll", "VmsyAll"
    )
  )
row_names <- paste("menhaden", ages)
row_names[length(row_names)] <- "menhaden 6+"

om_fmsy_group <- fmsy_data$Fbase[fmsy_data$Group %in% row_names] * compensation_msy$FMSY

# OM unfished biomass --------------------------------------------

# Only menhaden F = 0, other functional groups: F = stock assessment F from 1985 - 2017, then F = 0 for the rest of the years
biomass <- read_ewe_output(
  file_path = here::here("data", "ewe", "7ages_ecosim_om_unfished", "ecosim_forcing_pdsi_egg_amo1_F0_menhaden_1985"),
  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_f0_menhaden <- apply(biomass[[1]][, age_name], 1, sum)[time_id] * 1000000

biomass <- read_ewe_output(
  file_path = here::here("data", "ewe", "7ages_ecosim_om_unfished_base", "ecosim_base_run"),
  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_f0_menhaden_basecase <- apply(biomass[[1]][, age_name], 1, sum)[time_id] * 1000000

# Only menhaden and bass F = 0, other functional groups: F = stock assessment F from 1985 - 2017, then F = 0 for the rest of the years
biomass <- read_ewe_output(
  file_path = here::here("data", "ewe", "7ages_ecosim_om_unfished", "ecosim_forcing_pdsi_egg_amo1_F0_menhaden_bass_1985"),
  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_f0_menhaden_bass <- apply(biomass[[1]][, age_name], 1, sum)[time_id] * 1000000

# All functional groups F = 0: from 1985
biomass <- read_ewe_output(
  file_path = here::here("data", "ewe", "7ages_ecosim_om_unfished", "ecosim_forcing_pdsi_egg_amo1_F0_all_1985"),
  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_f0_all_1985 <- apply(biomass[[1]][, age_name], 1, sum)[time_id] * 1000000

x <- start_year:(start_year + length(time_id) - 1)
yrange <- range(biomass_f0_menhaden, biomass_f0_menhaden_bass, biomass_f0_all_1985)
plot(x, biomass_f0_menhaden,
  xlab = "Year", ylab = "Biomass",
  pch = 1, lty = 1, type = "l"
)
lines(x, biomass_f0_menhaden_bass,
  pch = 2, lty = 2
)
lines(x, biomass_f0_all_1985,
  pch = 3, lty = 3
)
legend("topright",
  c(
    "Menhaden F = 0",
    "Menhaden and bass F = 0",
    "All functional groups F = 0"
  ),
  pch = c(1, 2, 3),
  lty = c(1, 2, 3),
  bty = "n"
)

B0_F0_menhaden <- mean(tail(biomass_f0_menhaden, n = length(ages) * 2))
B0_F0_menhaden_basecase <- mean(tail(biomass_f0_menhaden_basecase, n = length(ages) * 2))
B0_F0_menhaden_bass <- mean(tail(biomass_f0_menhaden_bass, n = length(ages) * 2))
B0_F0_all <- mean(tail(biomass_f0_all_1985, n = length(ages) * 2))
B0_F0_mean <- mean(c(B0_F0_menhaden, B0_F0_menhaden_bass, B0_F0_all))

file_path <- ewe_output_path
set.seed(999)
# Load EwE biomass data

biomass <- read_ewe_output(
  file_path = file_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)[1:(length(model_year) + length(projection_year))]
biomass_ewe <- apply(biomass[[1]][, age_name], 1, sum) * 1000000

# source(here::here("Rscript", "load_indices.R"))
source(here::here("Rscript", "load_indices_alt_regression.R"))

output_dir <- here::here("data", "data_moderate", ewe_scenario_name, terminal_year)
if (!dir.exists(output_dir)) dir.create(output_dir, recursive = T)

# Set up JABBA base case assessment
# case 0
# BmsyK = 0.4 # shape = 1.2
# BmsyK = 0.45 # shape = 1.5
# BmsyK = 0.5 # shape = 2
# BmsyK <- stationary_msy$BMSY / B0_F0_menhaden
BmsyK = sum(compensation_msy$BMSY)*1000000/B0_F0_menhaden # shape = 3.23
# BmsyK = msy_mean["BMSY"]/B0_F0_mean # shape = 7.6
## Do not use effort data
ss_case0_input_noEffort <- generate_jabba(
  assessment_name = "case0_noEffort",
  output_dir = output_dir,
  sa_data = sa_data,
  BmsyK = BmsyK,
  model_year = model_year,
  projection_year = projection_year,
  effort_data = NULL
)


# k.prior (mu.k, cv.k)
if (add_environmental_effects == FALSE) {
  K.init <- B0_F0_menhaden_basecase / 1.45
} else {
  K.init <- B0_F0_menhaden
}

# K.init <- 8*max(sa_data$fishery$obs_total_catch_biomass$fleet1)
K.prior <- c(K.init, max(sa_data$fishery$obs_total_catch_biomass$fleet1))
log.K <- mean(log(K.prior))
sd.K <- abs(log.K - log(K.prior[1])) / 2
log.K <- mean(log(K.prior)) + 0.5 * sd.K^2 # Will be back-bias corrected in lognorm function
CV.K <- sqrt(exp(sd.K^2) - 1) #
# r.prior (mu.r, sd.r) # mu.r range: 0.2 - 0.8
start.r <- c(0.2, 0.8)
# start.r <- c(0.1, 0.4)
log.r <- mean(log(start.r))
sd.r <- abs(log.r - log(start.r[1])) / 2


# ss_case0_input_noEffort$jagsdata$K.pr <- c(k_init, 0.1)
# ss_case0_input_noEffort$jagsdata$r.pr <- c(0.2, 0.1)
# # ss_case0_input_noEffort$jagsdata$r.pr <- c(0.6, 0.1)
# psi_init <- biomass_ewe[1]/k_init
# ss_case0_input_noEffort$jagsdata$psi.pr <- c(psi_init, 0.1)

# k_init <- B0_F0_all
# k_init <- B0_F0_menhaden
ss_case0_input_noEffort$jagsdata$K.pr <- c(exp(log.K), CV.K) # 5208722, 2091
# ss_case0_input_noEffort$jagsdata$K.pr <- c(k_init, 0.5)
ss_case0_input_noEffort$jagsdata$r.pr <- c(exp(log.r), sd.r) # 0.4, 0.35
# psi_init <- biomass_ewe[1]/k_init
# psi_init <- 0.9
psi_init <- biomass_ewe[1] / exp(log.K)
ss_case0_input_noEffort$jagsdata$psi.pr <- c(psi_init, 0.25)

ss_case0_noEffort <- JABBA::fit_jabba(
  jbinput = ss_case0_input_noEffort,
  save.jabba = TRUE,
  save.all = TRUE,
  save.prj = TRUE,
  output.dir = output_dir,
  save.csvs = T
)

write.csv(ss_case0_noEffort$estimates, file.path(output_dir, "case0_noEffort_estimates.csv"))


## Use effort data
year_id <- seq(1, nrow(amo_unsmooth_lag1), by = 12)

if (add_fleet_dynamics == TRUE) {
  effort <- exp(menhadenE)
} else {
  effort <- menhadenE
}

ss_case0_effort_input <- generate_jabba(
  assessment_name = "case0",
  output_dir = output_dir,
  sa_data = sa_data,
  BmsyK = BmsyK,
  model_year = model_year,
  projection_year = projection_year,
  effort_data = effort[1:length(model_year)]
)

# k_init <- max(sa_data$fishery$obs_total_catch_biomass$fleet1) * 10
# k_init <- 3500000
# k_init <- B0_F0_menhaden
# ss_case0_effort_input$jagsdata$K.pr <- c(k_init, 0.1)
# ss_case0_effort_input$jagsdata$r.pr <- c(0.2, 0.1)
# # ss_case0_effort_input$jagsdata$r.pr <- c(0.6, 0.1)
# psi_init <- biomass_ewe[1]/k_init
# ss_case0_effort_input$jagsdata$psi.pr <- c(psi_init, 0.1)

# k_init <- B0_F0_all
# k_init <- B0_F0_menhaden
ss_case0_effort_input$jagsdata$K.pr <- c(exp(log.K), CV.K) # 5208722, 2091
# ss_case0_input_noEffort$jagsdata$K.pr <- c(k_init, 0.5)
ss_case0_effort_input$jagsdata$r.pr <- c(exp(log.r), sd.r) # 0.4, 0.35
# psi_init <- biomass_ewe[1]/k_init
psi_init <- biomass_ewe[1] / exp(log.K)
ss_case0_effort_input$jagsdata$psi.pr <- c(psi_init, 0.25)

set.seed(999)
ss_case0_effort <- JABBA::fit_jabba(
  jbinput = ss_case0_effort_input,
  save.jabba = TRUE,
  save.all = TRUE,
  save.prj = TRUE,
  output.dir = output_dir,
  save.csvs = T
)

write.csv(ss_case0_effort$estimates, file.path(output_dir, "case0_estimates.csv"))

ss_case0_input <- ss_case0_effort_input
ss_case0 <- ss_case0_effort

# Plot figures using JABBA
figure_dir <- file.path(figure_path)
if (!dir.exists(figure_dir)) dir.create(figure_dir)
jabba_figure_path <- file.path(figure_dir, terminal_year, "noEffort")
if (!dir.exists(jabba_figure_path)) dir.create(jabba_figure_path, recursive = T)
JABBA::jabba_plots(jabba = ss_case0_noEffort, output.dir = jabba_figure_path)
# plot with CIs (80% for projections)
# JABBA::jbplot_prj(ss_case0_noEffort, type = "BBmsy", output.dir = jabba_figure_path)
# JABBA::jbplot_prj(ss_case0_noEffort, type = "BB0", output.dir = jabba_figure_path)
# JABBA::jbplot_prj(ss_case0_noEffort, type = "FFmsy", output.dir = jabba_figure_path)

jabba_figure_path <- file.path(figure_dir, terminal_year, "Effort")
if (!dir.exists(jabba_figure_path)) dir.create(jabba_figure_path, recursive = T)
JABBA::jabba_plots(jabba = ss_case0_effort, output.dir = jabba_figure_path)
# plot with CIs (80% for projections)
# JABBA::jbplot_prj(ss_case0_effort, type = "BBmsy", output.dir = jabba_figure_path)
# JABBA::jbplot_prj(ss_case0_effort, type = "BB0", output.dir = jabba_figure_path)
# JABBA::jbplot_prj(ss_case0_effort, type = "FFmsy", output.dir = jabba_figure_path)

# Plot figures for cases 0.1 and 0.2 ----------------------------------------------------
fmsy0 <- rep(ss_case0_noEffort$posteriors$Hmsy, times = (length(projection_year) + 1) * length(ss_case0_input_noEffort$jagsdata$TAC[1, ]))
ss_case0_noEffort$prj_posterior$f <- ss_case0_noEffort$prj_posterior$harvest * fmsy0

bmsy0 <- rep(ss_case0_noEffort$posteriors$SBmsy, times = (length(projection_year) + 1) * length(ss_case0_input_noEffort$jagsdata$TAC[1, ]))
ss_case0_noEffort$prj_posterior$biomass <- ss_case0_noEffort$prj_posterior$stock * bmsy0

jpeg(filename = file.path(figure_path, paste0("Biomass", terminal_year, ewe_scenario_name, scenario_filename, ".jpeg")), width = 200, height = 120, units = "mm", res = 1800)
# plot biomass over time
par(mfrow = c(1, 2))

ylim <- range(
  biomass_ewe[time_id],
  ss_case0_noEffort$timeseries[, , "B"],
  ss_case0_noEffort$prj_posterior$biomass[(round(ss_case0_noEffort$prj_posterior$tac) %in% round(mean(tail(ss_case0_input_noEffort$data$catch$Total, n = 3))))],
  ss_case0$timeseries[, , "B"],
  ss_case0$prj_posterior$biomass[(round(ss_case0$prj_posterior$tac) %in% round(mean(tail(ss_case0_input$data$catch$Total, n = 3))))]
)

# Plot biomass over time based on no effort model
plot(c(model_year, projection_year),
  biomass_ewe[time_id],
  xlab = "Year", ylab = "Biomass (mt)",
  ylim = ylim,
  pch = 16
)

lines(model_year, ss_case0_noEffort$timeseries[, "mu", "B"])
lines(model_year, ss_case0_noEffort$timeseries[, "lci", "B"], lty = 2)
lines(model_year, ss_case0_noEffort$timeseries[, "uci", "B"], lty = 2)


# Plot biomass based on last three years of catch

for (i in 1:length(projection_year)) {
  # Projection TAC input used mean catch from the recent 3 years
  subdata_id <- (ss_case0_noEffort$prj_posterior$year == projection_year[i]) & (round(ss_case0_noEffort$prj_posterior$tac) %in% round(mean(tail(ss_case0_input_noEffort$data$catch$Total, n = 3))))

  boxplot(ss_case0_noEffort$prj_posterior$biomass[subdata_id],
    add = TRUE,
    at = projection_year[i],
    col = "gray",
    width = 1,
    outline = TRUE,
    axes = FALSE
  )
}
points(c(model_year, projection_year),
  biomass_ewe[time_id],
  pch = 16
)
legend("topleft", "No CPUE index", bty = "n")

# Plot biomass over time based model that is using effort data
fmsy0 <- rep(ss_case0$posteriors$Hmsy, times = (length(projection_year) + 1) * length(ss_case0_input$jagsdata$TAC[1, ]))
ss_case0$prj_posterior$f <- ss_case0$prj_posterior$harvest * fmsy0

bmsy0 <- rep(ss_case0$posteriors$SBmsy, times = (length(projection_year) + 1) * length(ss_case0_input$jagsdata$TAC[1, ]))
ss_case0$prj_posterior$biomass <- ss_case0$prj_posterior$stock * bmsy0

plot(c(model_year, projection_year),
  biomass_ewe[time_id],
  xlab = "Year", ylab = "Biomass (mt)",
  ylim = ylim,
  pch = 16
)

lines(model_year, ss_case0$timeseries[, "mu", "B"])
lines(model_year, ss_case0$timeseries[, "lci", "B"], lty = 2)
lines(model_year, ss_case0$timeseries[, "uci", "B"], lty = 2)

cor(biomass_ewe[time_id][1:length(model_year)], ss_case0$timeseries[, "mu", "B"])

# Plot biomass based on last three years of catch
for (i in 1:length(projection_year)) {
  subdata_id <- (ss_case0$prj_posterior$year == projection_year[i]) & (round(ss_case0$prj_posterior$tac) %in% round(mean(tail(ss_case0_input$data$catch$Total, n = 3))))

  boxplot(ss_case0$prj_posterior$biomass[subdata_id],
    add = TRUE,
    at = projection_year[i],
    col = "gray",
    width = 1,
    outline = TRUE,
    axes = FALSE
  )
}
points(c(model_year, projection_year),
  biomass_ewe[time_id],
  pch = 16
)
legend("topleft", "CPUE index included", bty = "n")
dev.off()

# Plot mortality over time

mortality <- read.csv(file = here::here("data", "ewe", "7ages_ecosim_om", paste0(ewe_scenario_name, "_fatage.csv")))

if (add_fleet_dynamics == TRUE) {
  mortality_ewe_apical <- apply(mortality[, age_name], 1, max)
  mortality_ewe_average <- apply(mortality[, age_name], 1, mean)
}

if (add_fleet_dynamics == FALSE) {
  mortality_ewe_apical <- apply(mortality[, paste0("menhaden", ages)], 1, max)
  mortality_ewe_average <- apply(mortality[, paste0("menhaden", ages)], 1, mean)
}


jpeg(filename = file.path(figure_path, paste0("fishing_mortality", terminal_year, ewe_scenario_name, scenario_filename, ".jpeg")), width = 200, height = 120, units = "mm", res = 1800)
par(mfrow = c(1, 1))
ylim <- range(
  mortality_ewe_apical,
  mortality_ewe_average,
  ss_case0$timeseries[, "mu", "F"]
)
plot(c(model_year, projection_year),
  # mortality_ewe[time_id],
  mortality_ewe_apical[1:(length(model_year) + projection_nyear)],
  type = "o",
  xlab = "Year", ylab = "F",
  pch = 1,
  col = 1,
  lty = 1,
  ylim = ylim
)
lines(c(model_year, projection_year),
  mortality_ewe_average[1:(length(model_year) + projection_nyear)],
  pch = 2, type = "o", lty = 2, col = 2
)
lines(model_year, ss_case0$timeseries[, "mu", "F"], lty = 3, col = 3)
legend("topleft",
  c("EwE Apical F", "EwE Average F", "JABBA F"),
  bty = "n",
  pch = c(1, 2, NA),
  lty = 1:3,
  col = 1:3
)

cor(mortality_ewe_average[1:length(model_year)], ss_case0$timeseries[, "mu", "F"])
dev.off()

Cases 1-5 are based on the settings from case 0

Case 1: Link Atlantic Multidecadal Oscillation Index with menhaden biomass estimates and adjust projections: AMO is an indicator of climate conditions and would affect recruitment variability of menhaden-like species
Case 2: Link Palmer drought severity index with menhaden biomass estimates and adjust projections: PSDI is a long-term indicator of drought conditions and it reflects river discharge and precipitation
Case 3: Link biomass of Striped bass from the EwE with menhaden biomass estimates and adjust projections because bass is a major predator
Case 4: Link fishing effort of menhaden with menhaden biomass estimates and adjust projections
Case 5: Link catch per unit effort of menhaden with menhaden biomass estimates and adjust projections
Linear regression models from case 1 - 5
True biomass of menhaden-like species as functions of AMO, PDSI, biomass of striped bass, fishing effort of menhaden, and menhaden CPUE

Status of indicators (SOI)

If stock-indicator relationship is positive, $SOI_{y} = \frac{I_{y}-I_{y}^{min}}{I_{y}^{max}-I_{y}^{min}}$
If stock-indicator relationship is negative, $SOI_{y} = 1-\frac{I_{y}-I_{y}^{min}}{I_{y}^{max}-I_{y}^{min}}$

where $I_{y}$ represents indicator value in year y. $I_{y}^{min}$ and $I_{y}^{max}$ represent the minimum and maximum values of $I$ from the time series.

Adjust projections

If $\frac{B2008_{i}}{BMSY} > 1$, $F^{'}_{i} = FMSY^{min} + SOI2008 \times (FMSY^{max}-FMSY^{min})$
If $\frac{B2008_{i}}{BMSY} \le 1$ and $\frac{B2008_{i}}{BMSY} > 0.5$, $F^{'}{i} = SOI2008 \times \frac{B2008{i}}{BMSY} \times FMSY_{i}$
If $\frac{B2008_{i}}{BMSY} \le 0.5$, $F^{'}_{i} = 0$

where $i$ represents individual iterations

# Projection: cases 1 - 5 ---------------------------

lm_slope <- data.frame(
  case = "PDSI+AMO1 Case",
  amo = NA,
  pdsi = NA,
  bassB = NA,
  menhadenC = NA,
  menhadenE = NA,
  menhadenCPUE = NA,
  bassCPUE = NA,
  herringCPUE = NA,
  menhadenV = NA
)

for (projection_year_id in 1:length(projection_year)) {
  menhaden_b <- ss_case0$timeseries[, "mu", "B"]

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

  index_year <- model_year


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

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

  y <- log(menhaden_b[biomass_lag_id])

    # Linear regression model -----------------------------------------
    # amo
    amo <- amo_unsmooth_lag1[year_id, ]
    amo_x <- amo$scaled_value[index_lag_id]
    lm_slope$amo <- regression_analysis(x = amo_x, y)[[1]]
    amo_model <- regression_analysis(x = amo_x, y)[[2]]

    # pdsi
    pdsi <- palmer_drought_severity_index[year_id, ]
    pdsi_x <- pdsi$scaled_value[index_lag_id]
    lm_slope$pdsi <- regression_analysis(x = pdsi_x, y)[[1]]
    pdsi_model <- regression_analysis(x = pdsi_x, y)[[2]]

    # bassB
    bassB <- bass_bio[bass_bio$Year %in% year_id, ]
    bassB_x <- bassB$bass_bio[index_lag_id]
    lm_slope$bassB <- regression_analysis(x = bassB_x, y)[[1]]
    bassB_model <- regression_analysis(x = bassB_x, y)[[2]]

    # menhadenC
    menhadenC <- sa_data$fishery$obs_total_catch_biomass$fleet1[1:length(year_id)]
    menhadenC_x <- menhadenC[index_lag_id]
    lm_slope$menhadenC <- regression_analysis(x = menhadenC_x, y)[[1]]
    menhadenC_model <- regression_analysis(x = menhadenC_x, y)[[2]]

    # menhadenE
    if (add_fleet_dynamics == TRUE) {
      menhadenEffort <- effort_all[, "F.menhaden"]
      # menhadenE <- log(menhadenEffort[year_id])
      menhadenE <- menhadenEffort[year_id]
    } else {
      menhadenE <- apply(fatage_all[year_id, grep("AtlanticMenhaden", colnames(fatage_all))] / menhadenCatchability[year_id, grep("menhaden", colnames(menhadenCatchability))], 1, sum)
    }
    menhadenE_x <- menhadenE[index_lag_id]
    lm_slope$menhadenE <- regression_analysis(x = menhadenE_x, y)[[1]]
    menhadenE_model <- regression_analysis(x = menhadenE_x, y)[[2]]

    # menhadenCPUE
    if (add_fleet_dynamics == TRUE) {
      menhadenCPUE <- menhadenC / menhadenEffort[year_id]
    } else {
      menhadenCPUE <- menhadenC / menhadenE
    }
    menhadenCPUE_x <- menhadenCPUE[index_lag_id]
    lm_slope$menhadenCPUE <- regression_analysis(x = menhadenCPUE_x, y)[[1]]
    menhadenCPUE_model <- regression_analysis(x = menhadenCPUE_x, y)[[2]]

    # bassCPUE
    bassCatch <- apply(as.data.frame(catch[[1]][, grep("StripedBass", names(catch[[1]]))]), 1, sum) * 1000000
    bassC <- bassCatch[year_id]
    if (add_fleet_dynamics == TRUE) {
      bassEffort <- effort_all[, "F.striped.bass"]
      bassE <- bassEffort[year_id]
    }
    if (add_fleet_dynamics == FALSE) {
      bassE <- apply(fatage_all[year_id, grep("StripedBass", colnames(fatage_all))] / bassCatchability[year_id, grep("striped.bass", colnames(bassCatchability))], 1, sum)
    }

  bassCPUE <- bassC / bassE
    bassCPUE_x <- bassCPUE[index_lag_id]
    lm_slope$bassCPUE <- regression_analysis(x = bassCPUE_x, y)[[1]]
    bassCPUE_model <- regression_analysis(x = bassCPUE_x, y)[[2]]

    # herringCPUE
    herringCatch <- apply(as.data.frame(catch[[1]][, grep("AtlanticHerring", names(catch[[1]]))]), 1, sum) * 1000000
    herringC <- herringCatch[year_id]
    if (add_fleet_dynamics == TRUE) {
      herringEffort <- effort_all[, "F.herring"]
      herringE <- herringEffort[year_id]
    }
    if (add_fleet_dynamics == FALSE) {
      herringE <- apply(fatage_all[year_id, grep("AtlanticHerring", colnames(fatage_all))] / herringCatchability[year_id, grep("Atlantic.herring", colnames(herringCatchability))], 1, sum)
    }

    herringCPUE <- herringC / herringE
    herringCPUE_x <- herringCPUE[index_lag_id]
    lm_slope$herringCPUE <- regression_analysis(x = herringCPUE_x, y)[[1]]
    herringCPUE_model <- regression_analysis(x = herringCPUE_x, y)[[2]]

    # menhadenV
    menhadenV <- menhadenValue[year_id]
    menhadenV_x <- menhadenV[index_lag_id]
    lm_slope$menhadenV <- regression_analysis(x = menhadenV_x, y)[[1]]
    menhadenV_model <- regression_analysis(x = menhadenV_x, y)[[2]]

    if (projection_year_id == 1) {
      lm_data_em <- rbind(
        data.frame(
          Menhaden_Biomass = y,
          Variable = "AMO",
          Index_Value = amo_x
        ),
        data.frame(
          Menhaden_Biomass = y,
          Variable = "PDSI",
          Index_Value = pdsi_x
        ),
        data.frame(
          Menhaden_Biomass = y,
          Variable = "Biomass of Striped bass",
          Index_Value = bassB_x
        ),
        data.frame(
          Menhaden_Biomass = y,
          Variable = "Menhaden Catch",
          Index_Value = menhadenC_x
        ),
        data.frame(
          Menhaden_Biomass = y,
          Variable = "Menhaden fishing effort",
          Index_Value = menhadenE_x
        ),
        data.frame(
          Menhaden_Biomass = y,
          Variable = "Menhaden CPUE",
          Index_Value = menhadenCPUE_x
        ),
        data.frame(
          Menhaden_Biomass = y,
          Variable = "Bass CPUE",
          Index_Value = bassCPUE_x
        ),
        data.frame(
          Menhaden_Biomass = y,
          Variable = "Herring CPUE",
          Index_Value = herringCPUE_x
        ),
        data.frame(
          Menhaden_Biomass = y,
          Variable = "Menhaden Ex-vessel Value",
          Index_Value = menhadenV_x
        )
      )
      lm_data_em$model <- "EM"
    }

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

    amo_soi <- calc_soi(
      # indicator_data = data_transformation(amo$scaled_value),
      indicator_data = amo$scaled_value,
      slope = coef(amo_model)["x_vec"]
    )

    pdsi_soi <- calc_soi(
      # indicator_data = data_transformation(pdsi$scaled_value),
      indicator_data = pdsi$scaled_value,
      slope = coef(pdsi_model)["x_vec"]
    )

    bassB_soi <- calc_soi(
      # indicator_data = data_transformation(bassB$bass_bio),
      indicator_data = bassB$bass_bio,
      slope = coef(bassB_model)["x_vec"]
    )

    menhadenC_soi <- calc_soi(
      # indicator_data = data_transformation(menhadenC),
      indicator_data = menhadenC,
      slope = coef(menhadenC_model)["x_vec"]
    )

    menhadenE_soi <- calc_soi(
      # indicator_data = data_transformation(menhadenE),
      indicator_data = menhadenE,
      slope = coef(menhadenE_model)["x_vec"]
    )

    menhadenCPUE_soi <- calc_soi(
      # indicator_data = data_transformation(menhadenCPUE),
      indicator_data = menhadenCPUE,
      slope = coef(menhadenCPUE_model)["x_vec"]
    )

    bassCPUE_soi <- calc_soi(
      # indicator_data = data_transformation(bassCPUE),
      indicator_data = bassCPUE,
      slope = coef(bassCPUE_model)["x_vec"]
    )

    herringCPUE_soi <- calc_soi(
      # indicator_data = data_transformation(herringCPUE),
      indicator_data = herringCPUE,
      slope = coef(herringCPUE_model)["x_vec"]
    )

    menhadenV_soi <- calc_soi(
      # indicator_data = data_transformation(menhadenV),
      indicator_data = menhadenV,
      slope = coef(menhadenV_model)["x_vec"]
    )

  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],
      menhadenC = scale(menhadenC)[, 1],
      menhadenE = scale(menhadenE)[, 1],
      menhadenCPUE = scale(menhadenCPUE)[, 1],
      bassCPUE = scale(bassCPUE)[, 1],
      herringCPUE = scale(herringCPUE)[, 1],
      menhadenV = scale(menhadenV)[, 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,
      bassB = bassB_soi,
      menhadenC = menhadenC_soi,
      menhadenE = menhadenE_soi,
      menhadenCPUE = menhadenCPUE_soi,
      bassCPUE = bassCPUE_soi,
      herringCPUE = herringCPUE_soi,
      menhadenV = menhadenV_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],
        menhadenC = scale(menhadenC)[, 1],
        menhadenE = scale(menhadenE)[, 1],
        menhadenCPUE = scale(menhadenCPUE)[, 1],
        bassCPUE = scale(bassCPUE)[, 1],
        herringCPUE = scale(herringCPUE)[, 1],
        menhadenV = scale(menhadenV)[, 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,
        bassB = bassB_soi,
        menhadenC = menhadenC_soi,
        menhadenE = menhadenE_soi,
        menhadenCPUE = menhadenCPUE_soi,
        bassCPUE = bassCPUE_soi,
        herringCPUE = herringCPUE_soi,
        menhadenV = menhadenV_soi
      )
    )
  }


  # Adjusted FMSY ----------------------------------------------------
  Bt_BMSY <- ss_case0$posteriors$BtoBmsy[, length(model_year)]

  amo_fmsy <- adjust_projection_jabba(
    FMSY = ss_case0$refpts_posterior$Fmsy,
    soi = tail(amo_soi, n = 1),
    Bt_BMSY = as.matrix(Bt_BMSY)
  )
  pdsi_fmsy <- adjust_projection_jabba(
    FMSY = ss_case0$refpts_posterior$Fmsy,
    soi = tail(pdsi_soi, n = 1),
    Bt_BMSY = as.matrix(Bt_BMSY)
  )
  bassB_fmsy <- adjust_projection_jabba(
    FMSY = ss_case0$refpts_posterior$Fmsy,
    soi = tail(bassB_soi, n = 1),
    Bt_BMSY = as.matrix(Bt_BMSY)
  )
  menhadenC_fmsy <- adjust_projection_jabba(
    FMSY = ss_case0$refpts_posterior$Fmsy,
    soi = tail(menhadenC_soi, n = 1),
    Bt_BMSY = as.matrix(Bt_BMSY)
  )
  menhadenE_fmsy <- adjust_projection_jabba(
    FMSY = ss_case0$refpts_posterior$Fmsy,
    soi = tail(menhadenE_soi, n = 1),
    Bt_BMSY = as.matrix(Bt_BMSY)
  )
  menhadenCPUE_fmsy <- adjust_projection_jabba(
    FMSY = ss_case0$refpts_posterior$Fmsy,
    soi = tail(menhadenCPUE_soi, n = 1),
    Bt_BMSY = as.matrix(Bt_BMSY)
  )
  bassCPUE_fmsy <- adjust_projection_jabba(
    FMSY = ss_case0$refpts_posterior$Fmsy,
    soi = tail(bassCPUE_soi, n = 1),
    Bt_BMSY = as.matrix(Bt_BMSY)
  )
  herringCPUE_fmsy <- adjust_projection_jabba(
    FMSY = ss_case0$refpts_posterior$Fmsy,
    soi = tail(herringCPUE_soi, n = 1),
    Bt_BMSY = as.matrix(Bt_BMSY)
  )
  menhadenV_fmsy <- adjust_projection_jabba(
    FMSY = ss_case0$refpts_posterior$Fmsy,
    soi = tail(menhadenV_soi, n = 1),
    Bt_BMSY = as.matrix(Bt_BMSY)
  )

  if (projection_year_id == 1) {
    fmsy_data <- data.frame(
      iter = 1:length(amo_fmsy),
      projection_year_id = projection_year[projection_year_id],
      JABBA = ss_case0$refpts_posterior$Fmsy,
      amo = amo_fmsy,
      pdsi = pdsi_fmsy,
      bassB = bassB_fmsy,
      menhadenC = menhadenC_fmsy,
      menhadenE = menhadenE_fmsy,
      menhadenCPUE = menhadenCPUE_fmsy,
      bassCPUE = bassCPUE_fmsy,
      herringCPUE = herringCPUE_fmsy,
      menhadenV = menhadenV_fmsy
    )
  } else {
    fmsy_data <- rbind(
      fmsy_data,
      data.frame(
        iter = 1:length(amo_fmsy),
        projection_year_id = projection_year[projection_year_id],
        JABBA = ss_case0$refpts_posterior$Fmsy,
        amo = amo_fmsy,
        pdsi = pdsi_fmsy,
        bassB = bassB_fmsy,
        menhadenC = menhadenC_fmsy,
        menhadenE = menhadenE_fmsy,
        menhadenCPUE = menhadenCPUE_fmsy,
        bassCPUE = bassCPUE_fmsy,
        herringCPUE = herringCPUE_fmsy,
        menhadenV = menhadenV_fmsy
      )
    )
  }
}

output <- list(
  lm_data_em = lm_data_em,
  lm_slope = lm_slope,
  soi_data = soi_data,
  fmsy_data = fmsy_data,
  Bt_BMSY = as.matrix(Bt_BMSY)
)

output$lm_data_om <- lm_data_om

# Linear regression figure
lm_data <- rbind(lm_data_om, lm_data_em)
lm_data$Variable <- factor(lm_data$Variable, levels = c("AMO", "PDSI", "Biomass of Striped bass", "Menhaden Catch", "Menhaden fishing effort", "Menhaden CPUE", "Bass CPUE", "Herring CPUE", "Menhaden Ex-vessel Value"))

ggplot(lm_data, aes(x = Index_Value, y = Menhaden_Biomass, color = model)) +
  geom_point() +
  geom_smooth(method = lm) +
  facet_wrap(~Variable, scales = "free") +
  # facet_wrap(~Variable, scales = "free_x") +
  labs(
    title = "Linear regression analysis",
    x = "Indicator Value",
    y = "Menhaden biomass (mt)"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "_linear_regression.jpeg")))

soi_data_melt <- reshape2::melt(
  soi_data,
  id = c("year", "projection_year_id")
)
soi_data_melt$projection_year_id <- rep(rep(projection_year, each = length(model_year)), times = ncol(soi_data) - 2)

output$soi_data_melt <- soi_data_melt

# SOI figure
ggplot(soi_data_melt[soi_data_melt$projection_year_id == projection_year[1], ], aes(x = year, y = value)) +
  geom_line() +
  facet_wrap(~variable, scales = "free_y") +
  labs(
    title = "Status of Indicators",
    x = "Year",
    y = "Status of Indicators"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "_soi.jpeg")))


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)

output$fmsy_data_melt <- fmsy_data_melt

ggplot(fmsy_data_melt, aes(x = iter, y = value, color = projection_year_id)) +
  geom_line() +
  facet_wrap(~variable) +
  labs(
    title = "FMSY",
    x = "Iteration",
    y = "FMSY"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "_fmsy_iter.jpeg")))

# FMSY figure
ggplot(fmsy_data_melt, aes(x = variable, y = value, color = projection_year_id)) +
  geom_boxplot() +
  labs(
    title = "FMSY",
    x = "",
    y = "FMSY"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "_FMSY.jpeg")))

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

output$fmsy_data_melt_median <- fmsy_data_melt_median

ggplot(fmsy_data_melt_median, aes(x = projection_year_id, y = value, color = variable)) +
  geom_line() +
  labs(
    title = "Median FMSY from JABBA and adjusted FMSY",
    x = "",
    y = "FMSY"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "_FMSY_median.jpeg")))

projection_f <- rbind(
  fmsy_data_melt_median,
  data.frame(
    projection_year_id = projection_year,
    variable = "om",
    value = max(compensation_msy$FMSY)
  )
)

fishery_sel <- IP4EBFM::logistic(
  pattern = "double_logistic",
  x = sa_data$biodata$ages,
  slope_asc = 3.1,
  location_asc = 1.8,
  slope_desc = 0.88,
  location_desc = 0.01
)

if (add_fleet_dynamics == TRUE) fishery_sel <- c(0.078, 0.17, 0.92, 0.54, 1, 0.51, 0.78) # From SS3

for (i in 1:nrow(projection_f)) {
  projection_f[i, paste0("fmsy_", sa_data$biodata$ages)] <- projection_f$value[i] * fishery_sel
  projection_f[i, paste0("fmsy1.25_", sa_data$biodata$ages)] <- projection_f$value[i] * fishery_sel * 1.25
  projection_f[i, paste0("fmsy0.75_", sa_data$biodata$ages)] <- projection_f$value[i] * fishery_sel * 0.75

  # if (projection_f[i, "variable"] == "om" & !is.null(om_fmsy_group)) {
  #   projection_f[i, paste0("fmsy_", sa_data$biodata$ages)] <- om_fmsy_group
  #   projection_f[i, paste0("fmsy1.25_", sa_data$biodata$ages)] <- om_fmsy_group*1.25
  #   projection_f[i, paste0("fmsy0.75_", sa_data$biodata$ages)] <- om_fmsy_group*0.75
  #
  # } else {
  #   projection_f[i, paste0("fmsy_", sa_data$biodata$ages)] <- projection_f$value[i] * fishery_sel
  #   projection_f[i, paste0("fmsy1.25_", sa_data$biodata$ages)] <- projection_f$value[i] * fishery_sel*1.25
  #   projection_f[i, paste0("fmsy0.75_", sa_data$biodata$ages)] <- projection_f$value[i] * fishery_sel*0.75
  # }
}
projection_f[which(projection_f$variable == "om"), 3:ncol(projection_f)] <- data.frame(t(c(max(om_fmsy_group), rep(om_fmsy_group, times=3))))
write.csv(fmsy_data_melt, here::here("data", "data_moderate", ewe_scenario_name, paste0("fmsy_data_", terminal_year, scenario_filename, ".csv")))
write.csv(projection_f, here::here("data", "data_moderate", ewe_scenario_name, paste0("fmsy_median_data_", terminal_year, scenario_filename, ".csv")))

Slope values from linear regression models

knitr::kable(lm_slope)
write.csv(output$lm_slope, here::here("data", "data_moderate", ewe_scenario_name, paste0("jabba_em_lm_slope_", terminal_year, scenario_filename, ".csv")))

Adjust catch input and project biomass based on FMSY*B

Bprojection_startyr <- ss_case0$posteriors$BtoBmsy[, (length(model_year) + 1)] * ss_case0$posteriors$SBmsy[, 1]

variable <- unique(fmsy_data_melt$variable)

Bt <- rep(Bprojection_startyr, times = length(variable))

tac_data_melt <- fmsy_data_melt
names(tac_data_melt)[names(tac_data_melt) == "value"] <- "fmsy"
tac_data_melt$tac <- adjust_projection_catcheco_jabba(fmsy_melted = fmsy_data_melt, Bt = Bt)
tac_data_melt_median <- aggregate(tac ~ variable, data = tac_data_melt, median)

output$tac_data_melt <- tac_data_melt
output$tac_data_melt_median <- tac_data_melt_median
save(output, file = here::here("data", "data_moderate", ewe_scenario_name, paste0("jabba_output_", terminal_year, scenario_filename, ".RData")))

ss_case0_input$jagsdata$nTAC <- length(variable)
ss_case0_input$jagsdata$TAC <- matrix(rep(tac_data_melt_median$tac, times = 5),
  ncol = length(variable), byrow = T
)
ss_case0_input$settings$TAC.implementation <- projection_year[1]
ss_case0 <- JABBA::fit_jabba(
  jbinput = ss_case0_input,
  save.jabba = TRUE,
  save.all = TRUE,
  save.prj = TRUE,
  output.dir = output_dir,
  save.csvs = T
)

biomass_projection <- ss_case0$prj_posterior
biomass_projection <- biomass_projection[(biomass_projection$year %in% projection_year), ]
Bmsy <- rep(ss_case0$posteriors$SBmsy[, 1],
  times = length(unique(biomass_projection$tac)) * length(unique(biomass_projection$year))
)
biomass_projection$biomass <- biomass_projection$stock * Bmsy
biomass_projection$year <- as.factor(biomass_projection$year)
biomass_projection$tac <- as.factor(biomass_projection$tac)

biomass_projection <- merge(biomass_projection, tac_data_melt_median, by = "tac")

# jpeg(filename = file.path(figure_path, paste0("biomass_projection", terminal_year, scenario_filename, ".jpeg")), width = 200, height = 100, units = "mm", res = 1800)
ggplot(biomass_projection, aes(x = year, y = biomass, color = year)) +
  geom_boxplot() +
  facet_wrap(~variable) +
  labs(
    title = " Estimated biomass based on catch=FMSY*B",
    x = "",
    y = "Biomass"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "_jabba_biomass_projection.jpeg")))

om_K <- B0_F0_menhaden
om_bmsy <- sum(compensation_msy$BMSY)
om_fmsy <- max(compensation_msy$FMSY)

cases <- c(
  "om", "jabba",
  "amo", "pdsi",
  "bassb",
  "menhadenc", "menhadene", "menhadencpue",
  "basscpue", "herringcpue",
  "menhadenv"
)

evaluation_table <- data.frame(
  Average_Catch = NA,
  B_K = NA,
  B_BMSY = NA,
  Average_Biomass = NA,
  Bonanza_Period = NA,
  Collapse_Period = NA
)

for (i in seq_along(cases)) {
  if (add_fleet_dynamics == FALSE) {
    file_path <- here::here("data", "ewe", "7ages_ecosim_om_projections", paste0(ewe_scenario_name, "_data_moderate_", terminal_year, scenario_filename, "_", cases[i]))
  }

  if (add_fleet_dynamics == TRUE) {
    file_path <- here::here("data", "ewe", "7ages_ecosim_om_projections_fleet_dynamics", paste0(ewe_scenario_name, "_annualF_data_moderate_", terminal_year, scenario_filename, "_", cases[i]))
  }

  # Load EwE biomass data
  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 <- read_ewe_output(
    file_path = file_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)[1:(length(model_year) + projection_nyear)]
  biomass_ewe <- apply(biomass[[1]][, age_name], 1, sum) * 1000000

  catch <- read_ewe_output(
    file_path = file_path,
    file_names = "catch_monthly.csv",
    skip_nrows = 8,
    plot = FALSE,
    figure_titles = NULL,
    functional_groups = functional_groups,
    figure_colors = NULL
  )
  catch_ewe <- apply(catch[[1]][, age_name], 1, sum) * 1000000

  evaluation_table[i, ] <- c(
    # mean catch
    mean(tail(catch_ewe[time_id], n = length(projection_year))),
    # Bend/K
    tail(biomass_ewe[time_id], n = 1) / om_K,
    # Bend/BMSY
    tail(biomass_ewe[time_id], n = 1) / om_bmsy,
    # mean biomass
    mean(tail(biomass_ewe[time_id], n = length(projection_year))),
    # Bonanza length
    sum((tail(biomass_ewe[time_id], n = length(projection_year)) / om_K) > 0.8),
    # Minmize collapse length
    sum((tail(biomass_ewe[time_id], n = length(projection_year)) / om_K) < 0.2)
  )
}

case_names <- levels(biomass_projection$variable)
for (i in 1:length(levels(biomass_projection$variable))) {
  projected_biomass <- aggregate(biomass ~ year + tac + variable, data = biomass_projection, median)
  bk <- aggregate(bk ~ year + tac + variable, data = biomass_projection, median)

  evaluation_table <- rbind(
    evaluation_table,
    c(
      # mean catch
      as.numeric(as.character(unique(biomass_projection$tac[biomass_projection$variable == case_names[i]]))),
      # median(biomass_projection$harvest[biomass_projection$year == tail(projection_year, n = 1) &
      #   biomass_projection$tac == unique(biomass_projection$tac[biomass_projection$variable == case_names[i]])]),
      # Bend/K
      median(biomass_projection$bk[biomass_projection$year == tail(projection_year, n = 1) &
        biomass_projection$tac == unique(biomass_projection$tac[biomass_projection$variable == case_names[i]])]),
      # Bend/BMSY
      median(biomass_projection$stock[biomass_projection$year == tail(projection_year, n = 1) &
        biomass_projection$tac == unique(biomass_projection$tac[biomass_projection$variable == case_names[i]])]),
      # Mean biomass
      mean(projected_biomass$biomass[projected_biomass$tac == unique(biomass_projection$tac[biomass_projection$variable == case_names[i]])]),
      # Bonanza length
      sum(bk$bk[bk$tac == unique(bk$tac[bk$variable == case_names[i]])] > 0.8),
      # Minimize collapse length
      sum(bk$bk[bk$tac == unique(bk$tac[bk$variable == case_names[i]])] < 0.2)
    )
  )
}

evaluation_table <- rbind(
  apply(evaluation_table, 2, max),
  rep(0, ncol(evaluation_table)),
  evaluation_table
)

row.names(evaluation_table) <- c(
  "Max", "Min",
  "OM+OM FMSY",
  "OM+DBSRA FMSY",
  "OM+AMO Fadj",
  "OM+PDSI Fadj",
  "OM+Bass Biomass Fadj",
  "OM+Menhaden Catch Fadj",
  "OM+Menhaden Effort Fadj",
  "OM+Menhaden CPUE Fadj",
  "OM+Bass CPUE Fadj",
  "OM+Herring CPUE Fadj",
  "OM+Menhaden Value Fadj",
  "DBSRA EM+DBSRA FMSY",
  "DBSRA EM+AMO Fadj",
  "DBSRA EM+PDSI Fadj",
  "DBSRA EM+Bass Biomass Fadj",
  "DBSRA EM+Menhaden Catch Fadj",
  "DBSRA EM+Menhaden Effort Fadj",
  "DBSRA EM+Menhaden CPUE Fadj",
  "DBSRA EM+Bass CPUE Fadj",
  "DBSRA EM+Herring CPUE Fadj",
  "DBSRA EM+Menhaden Value Fadj"
)

evaluation_table["Max", c("Bonanza_Period", "Collapse_Period")] <- c(length(projection_year), length(projection_year))
evaluation_table[c("Max", "Min"), c("Collapse_Period")] <- rev(evaluation_table[c("Max", "Min"), c("Collapse_Period")])

write.csv(evaluation_table, here::here("data", "data_moderate", ewe_scenario_name, paste0("evaluation_table_", terminal_year, scenario_filename, ".csv")))


case_num <- 10
jpeg(filename = file.path(figure_path, paste0(ewe_scenario_name, "_radar_mediumF_", terminal_year, scenario_filename, ".jpeg")), width = 200, height = 200, units = "mm", res = 1800)
par(mfrow = c(4, 3), mar = c(0, 0, 0, 0))
# colors <- c("gray", brewer.pal(12, "Paired"))

colors <- Polychrome::glasbey.colors(24)
colors <- colors[2:length(colors)]
# colors <- brewer.pal(12,"Paired")

variable_label <- c("Mean Catch", paste0("B", tail(projection_year, n = 1), "/K"), paste0("B", tail(projection_year, n = 1), "/BMSY"), "Mean Biomass", "Bonanza \n Length", "Minimize \n Collapse \n Length")

for (i in 1:case_num) {
  colors_border <- colors[c(1:2, 2 + i, case_num + 3, case_num + 3 + i)]
  radarchart(evaluation_table[c(1:4, 4 + i, case_num + 5, case_num + 5 + i), ],
    seg = 5,
    pcol = colors_border, plwd = 2, plty = 1,
    axistype = 0, axislabcol = "black",
    paxislabels = round(evaluation_table["Max", ], 2), palcex = 1,
    vlabels = variable_label, vlcex = 0.8
  )
  # colors_border <- colors[c(1, 1+i, case_num+2, case_num+2+i)]
  # radarchart(evaluation_table[c(1:3, 3+i, case_num+4, case_num+4+i), ],
  #            seg=5,
  #            pcol=colors_border, plwd=2, plty=1,
  #            axistype = 0, axislabcol="black",
  #            paxislabels = round(evaluation_table["Max",],2), palcex = 1,
  #            vlabels = variable_label, vlcex = 0.8)
}
colors_border <- colors[c(1:2, (case_num + 3):length(colors))]
radarchart(evaluation_table[c(1:4, (case_num + 5):nrow(evaluation_table)), ],
  seg = 5,
  pcol = colors_border, plwd = 2, plty = 1,
  axistype = 2, axislabcol = "black",
  paxislabels = round(evaluation_table["Max", ], 2), palcex = 1,
  vlabels = variable_label, vlcex = 0.8
)
colors_border <- colors[1:length(colors)]
plot(1, type = "n", axes = FALSE, xlab = "", ylab = "")
legend("center", legend = rownames(evaluation_table[-c(1, 2), ]), bty = "n", pch = 20, col = colors_border, text.col = "black", cex = 0.7, pt.cex = 1)

format(evaluation_table, digits = 3)
dev.off()

Bai-Li-NOAA/IFA4EBFM documentation built on July 1, 2024, 4:53 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 Projections for Ecosystem-Based Fisheries Management

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

Case 0: stock assessment base run (terminal year = 2008)

Cases 1-5 are based on the settings from case 0

Status of indicators (SOI)

Adjust projections

Adjust catch input and project biomass based on FMSY*B

R Package Documentation

Browse R Packages

We want your feedback!

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

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

Case 0: stock assessment base run (terminal year = 2008)

Cases 1-5 are based on the settings from case 0

Status of indicators (SOI)

Adjust projections

Adjust catch input and project biomass based on FMSY*B

R Package Documentation

Browse R Packages

We want your feedback!

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