In Bai-Li-NOAA/IFA4EBFM: Interdisciplinary Forecasting Approach 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

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_forcing_pdsi_egg_amo1")
} else {
  msy_file_path <- here::here("data", "ewe", "7ages_ecosim_om", "msy_base_run")
}


compensation_msy <- read_ewe_reference_points(
  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,
  plot = TRUE
)
# MSY, FMSY, BMSY: 0.5913518, 2.939998, 1.619142
# MSY, FMSY, BMSY: 591352, 2.9, 1619142

stationary_msy <- read_ewe_reference_points(
  file_path = msy_file_path,
  file_names = c(
    "F.menhaden_Catch_StationarySystem.csv",
    "F.menhaden_B_StationarySystem.csv"
  ),
  reference_points_scenario = "stationary",
  functional_groups = functional_groups,
  key_functional_group = "AtlanticMenhaden",
  ages = ages,
  plot = TRUE
)
# MSY, FMSY, BMSY: 0.5784879, 2.899998, 1.604036
# MSY, FMSY, BMSY: 578488, 2.9, 1604036

msy_mean <- sapply(1:length(compensation_msy), function(x) {
  mean(c(compensation_msy[[x]], stationary_msy[[x]]))
})
names(msy_mean) <- names(compensation_msy)

temp <- scan(file.path(msy_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$Fmsy[fmsy_data$Group %in% row_names]

# 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"))

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 = compensation_msy$BMSY/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)

  shapiro.test(menhaden_b[biomass_lag_id]) # < 0.05, data is not normally distributed
  shapiro.test(log(menhaden_b[biomass_lag_id])) # < 0.05, data is not normally distributed
  log_menhaden_b <- log(menhaden_b[biomass_lag_id])

  amo <- amo_unsmooth_lag1[year_id, ]
  amo_lm <- lm(log_menhaden_b ~ amo$scaled_value[index_lag_id])
  amo_fit <- fitted(amo_lm)
  amo_shapiro <- shapiro.test(summary(amo_lm)$residuals)
  lm_slope$amo <- paste0(
    ifelse(round(summary(amo_lm)$coefficients[2, 1], digits = 2) == 0,
      summary(amo_lm)$coefficients[2, 1],
      round(summary(amo_lm)$coefficients[2, 1], digits = 2)
    ),
    if (summary(amo_lm)$coefficients[2, 4] <= 0.05) {
      "*"
    },
    if (summary(amo_lm)$r.squared >= 0.6) {
      "^"
    },
    if (amo_shapiro$p.value > 0.05) {
      "~"
    }
  )

  pdsi <- palmer_drought_severity_index[year_id, ]
  pdsi_lm <- lm(log_menhaden_b ~ pdsi$scaled_value[index_lag_id])
  pdsi_fit <- fitted(pdsi_lm)
  pdsi_shapiro <- shapiro.test(summary(pdsi_lm)$residuals)
lm_slope$pdsi <- paste0(
  ifelse(round(summary(pdsi_lm)$coefficients[2, 1], digits = 2) == 0,
         summary(pdsi_lm)$coefficients[2, 1],
         round(summary(pdsi_lm)$coefficients[2, 1], digits = 2)
  ),
  if (summary(pdsi_lm)$coefficients[2, 4] <= 0.05) {
    "*"
  },
  if (summary(pdsi_lm)$r.squared >= 0.6) {
    "^"
  },
  if (pdsi_shapiro$p.value > 0.05) {
    "~"
  }
)

  # if (add_fleet_dynamics == TRUE){
  #   bassB <- log(bass_bio[bass_bio$Year %in% year_id, ])
  # } else {
  bassB <- bass_bio[bass_bio$Year %in% year_id, ]
  # }
  bassB_lm <- lm(log_menhaden_b ~ bassB$bass_bio[index_lag_id])
  bassB_fit <- fitted(bassB_lm)
  bassB_shapiro <- shapiro.test(summary(bassB_lm)$residuals)
lm_slope$bassB <- paste0(
  ifelse(round(summary(bassB_lm)$coefficients[2, 1], digits = 2) == 0,
         summary(bassB_lm)$coefficients[2, 1],
         round(summary(bassB_lm)$coefficients[2, 1], digits = 2)
  ),
  if (summary(bassB_lm)$coefficients[2, 4] <= 0.05) {
    "*"
  },
  if (summary(bassB_lm)$r.squared >= 0.6) {
    "^"
  },
  if (bassB_shapiro$p.value > 0.05) {
    "~"
  }
)

  # menhadenCatch <- apply(as.data.frame(catch[[1]][, grep("AtlanticMenhaden", names(catch[[1]]))]), 1, sum) *1000000

  # if (add_fleet_dynamics == TRUE){
  #   menhadenC <- log(menhadenCatch[year_id])
  # } else {
  # menhadenC <- menhadenCatch[year_id]
  # }
  menhadenC <- sa_data$fishery$obs_total_catch_biomass$fleet1[1:length(year_id)]
  shapiro.test(menhadenC) # < 0.05, data is not normally distributed
  shapiro.test(log(menhadenC)) # < 0.05, data is not normally distributed
  # menhadenC_lm <- lm(log_menhaden_b ~ menhadenC[index_lag_id])
  menhadenC_lm <- lm(log_menhaden_b ~ log(menhadenC[index_lag_id]))
  menhadenC_fit <- fitted(menhadenC_lm)
  bassB_shapiro <- shapiro.test(summary(bassB_lm)$residuals)
lm_slope$bassB <- paste0(
  ifelse(round(summary(bassB_lm)$coefficients[2, 1], digits = 2) == 0,
         summary(bassB_lm)$coefficients[2, 1],
         round(summary(bassB_lm)$coefficients[2, 1], digits = 2)
  ),
  if (summary(bassB_lm)$coefficients[2, 4] <= 0.05) {
    "*"
  },
  if (summary(bassB_lm)$r.squared >= 0.6) {
    "^"
  },
  if (bassB_shapiro$p.value > 0.05) {
    "~"
  }
)

  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_lm <- lm(log_menhaden_b ~ log(menhadenE[index_lag_id]))
  menhadenE_fit <- fitted(menhadenE_lm)
  menhadenE_shapiro <- shapiro.test(summary(menhadenE_lm)$residuals)
lm_slope$menhadenE <- paste0(
  ifelse(round(summary(menhadenE_lm)$coefficients[2, 1], digits = 2) == 0,
         summary(menhadenE_lm)$coefficients[2, 1],
         round(summary(menhadenE_lm)$coefficients[2, 1], digits = 2)
  ),
  if (summary(menhadenE_lm)$coefficients[2, 4] <= 0.05) {
    "*"
  },
  if (summary(menhadenE_lm)$r.squared >= 0.6) {
    "^"
  },
  if (menhadenE_shapiro$p.value > 0.05) {
    "~"
  }
)

  if (add_fleet_dynamics == TRUE) {
    menhadenCPUE <- menhadenC / menhadenEffort[year_id]
  } else {
    menhadenCPUE <- menhadenC / menhadenE
  }
  shapiro.test(menhadenCPUE) # < 0.05, data is not normally distributed
  shapiro.test(log(menhadenCPUE)) # = 0.5
  menhadenCPUE_lm <- lm(log_menhaden_b ~ log(menhadenCPUE[index_lag_id]))
  menhadenCPUE_fit <- fitted(menhadenCPUE_lm)
  menhadenCPUE_shapiro <- shapiro.test(summary(menhadenCPUE_lm)$residuals)
lm_slope$menhadenCPUE <- paste0(
  ifelse(round(summary(menhadenCPUE_lm)$coefficients[2, 1], digits = 2) == 0,
         summary(menhadenCPUE_lm)$coefficients[2, 1],
         round(summary(menhadenCPUE_lm)$coefficients[2, 1], digits = 2)
  ),
  if (summary(menhadenCPUE_lm)$coefficients[2, 4] <= 0.05) {
    "*"
  },
  if (summary(menhadenCPUE_lm)$r.squared >= 0.6) {
    "^"
  },
  if (menhadenCPUE_shapiro$p.value > 0.05) {
    "~"
  }
)

  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_lm <- lm(log_menhaden_b ~ bassCPUE[index_lag_id])
  bassCPUE_fit <- fitted(bassCPUE_lm)
  bassCPUE_shapiro <- shapiro.test(summary(bassCPUE_lm)$residuals)
lm_slope$bassCPUE <- paste0(
  ifelse(round(summary(bassCPUE_lm)$coefficients[2, 1], digits = 2) == 0,
         summary(bassCPUE_lm)$coefficients[2, 1],
         round(summary(bassCPUE_lm)$coefficients[2, 1], digits = 2)
  ),
  if (summary(bassCPUE_lm)$coefficients[2, 4] <= 0.05) {
    "*"
  },
  if (summary(bassCPUE_lm)$r.squared >= 0.6) {
    "^"
  },
  if (bassCPUE_shapiro$p.value > 0.05) {
    "~"
  }
)

  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_lm <- lm(log_menhaden_b ~ herringCPUE[index_lag_id])
  herringCPUE_fit <- fitted(herringCPUE_lm)
  herringCPUE_shapiro <- shapiro.test(summary(herringCPUE_lm)$residuals)
lm_slope$herringCPUE <- paste0(
  ifelse(round(summary(herringCPUE_lm)$coefficients[2, 1], digits = 2) == 0,
         summary(herringCPUE_lm)$coefficients[2, 1],
         round(summary(herringCPUE_lm)$coefficients[2, 1], digits = 2)
  ),
  if (summary(herringCPUE_lm)$coefficients[2, 4] <= 0.05) {
    "*"
  },
  if (summary(herringCPUE_lm)$r.squared >= 0.6) {
    "^"
  },
  if (herringCPUE_shapiro$p.value > 0.05) {
    "~"
  }
)

  # menhadenValue = value_all[, "F.menhaden"]
  # if (add_fleet_dynamics == TRUE){
  #   menhadenV <- log(menhadenValue[year_id])
  # } else {
  menhadenV <- menhadenValue[year_id]
  # }
  menhadenV_lm <- lm(log_menhaden_b ~ log(menhadenV[index_lag_id]))
  menhadenV_fit <- fitted(menhadenV_lm)
  menhadenV_shapiro <- shapiro.test(summary(menhadenV_lm)$residuals)
lm_slope$menhadenV <- paste0(
  ifelse(round(summary(menhadenV_lm)$coefficients[2, 1], digits = 2) == 0,
         summary(menhadenV_lm)$coefficients[2, 1],
         round(summary(menhadenV_lm)$coefficients[2, 1], digits = 2)
  ),
  if (summary(menhadenV_lm)$coefficients[2, 4] <= 0.05) {
    "*"
  },
  if (summary(menhadenV_lm)$r.squared >= 0.6) {
    "^"
  },
  if (menhadenV_shapiro$p.value > 0.05) {
    "~"
  }
)

  if (projection_year_id == 1) {
    lm_data_em <- rbind(
      data.frame(
        Menhaden_Biomass = log_menhaden_b,
        Variable = "AMO",
        Index_Value = amo$scaled_value[index_lag_id]
      ),
      data.frame(
        Menhaden_Biomass = log_menhaden_b,
        Variable = "PDSI",
        Index_Value = pdsi$scaled_value[index_lag_id]
      ),
      data.frame(
        Menhaden_Biomass = log_menhaden_b,
        Variable = "Biomass of Striped bass",
        Index_Value = bassB$bass_bio[index_lag_id]
      ),
      data.frame(
        Menhaden_Biomass = log_menhaden_b,
        Variable = "Menhaden Catch",
        Index_Value = log(menhadenC[index_lag_id])
      ),
      data.frame(
        Menhaden_Biomass = log_menhaden_b,
        Variable = "Menhaden fishing effort",
        Index_Value = log(menhadenE[index_lag_id])
      ),
      data.frame(
        Menhaden_Biomass = log_menhaden_b,
        Variable = "Menhaden CPUE",
        Index_Value = log(menhadenCPUE[index_lag_id])
      ),
      data.frame(
        Menhaden_Biomass = log_menhaden_b,
        Variable = "Bass CPUE",
        Index_Value = bassCPUE[index_lag_id]
      ),
      data.frame(
        Menhaden_Biomass = log_menhaden_b,
        Variable = "Herring CPUE",
        Index_Value = herringCPUE[index_lag_id]
      ),
      data.frame(
        Menhaden_Biomass = log_menhaden_b,
        Variable = "Menhaden Ex-vessel Value",
        Index_Value = log(menhadenV[index_lag_id])
      )
    )
    lm_data_em$model <- "EM"
  }

  # 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"
  #   )
  #   abline(amo_lm)
  #
  #   plot(pdsi$scaled_value[index_lag_id], menhaden_b[biomass_lag_id],
  #        xlab = "PDSI values",
  #        ylab = "Biomass of menhaden-like species"
  #   )
  #   abline(pdsi_lm)
  #
  #   plot(bassB$bass_bio[index_lag_id], menhaden_b[biomass_lag_id],
  #        xlab = "Biomass of Striped bass",
  #        ylab = "Biomass of menhaden-like species"
  #   )
  #   abline(bassB_lm)
  #
  #   plot(effort[index_lag_id], menhaden_b[biomass_lag_id],
  #        xlab = "Menhaden fishing effort",
  #        ylab = "Biomass of menhaden-like species"
  #   )
  #   abline(effort_lm)
  #
  # }


  # 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]
  )

  menhadenC_soi <- calc_soi(
    indicator_data = menhadenC,
    slope = coef(menhadenC_lm)[2]
  )

  menhadenE_soi <- calc_soi(
    indicator_data = menhadenE,
    slope = coef(menhadenE_lm)[2]
  )

  menhadenCPUE_soi <- calc_soi(
    indicator_data = menhadenCPUE,
    slope = coef(menhadenCPUE_lm)[2]
  )

  bassCPUE_soi <- calc_soi(
    indicator_data = bassCPUE,
    slope = coef(bassCPUE_lm)[2]
  )

  herringCPUE_soi <- calc_soi(
    indicator_data = herringCPUE,
    slope = coef(herringCPUE_lm)[2]
  )

  menhadenV_soi <- calc_soi(
    indicator_data = menhadenV,
    slope = coef(menhadenV_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],
      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 = compensation_msy$FMSY
  )
)

fishery_sel <- IFA4EBFM::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
  # }
}

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 <- compensation_msy$BMSY
om_fmsy <- 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 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

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 Forecasting Approach for Ecosystem-Based Fisheries Management

In Bai-Li-NOAA/IFA4EBFM: Interdisciplinary Forecasting Approach 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 Forecasting Approach for Ecosystem-Based Fisheries Management