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

Package index
Search the Bai-Li-NOAA/IFA4EBFM package
Vignettes
Functions
37
Source code
31
Man pages
18
Home
/
GitHub
/
Bai-Li-NOAA/IFA4EBFM
/

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

Case 0: stock assessment base run

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

# 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_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_poor", ewe_scenario_name)
if (!dir.exists(output_dir)) dir.create(output_dir)

# Scenario A: default stock assessment run ----------------------------
## DB-SRA is a method designed for determining a catch limit
## and management reference points for data-limited fisheries
## where catches are known from the beginning of exploitation.
## However, the catch data from the menhaden-like species case is not
## from the beginning of the exploitation.

# Populate the input data object

# Create a blank DLM object
ss_case0 <- new("Data")
# Change default area from 2 to 1
ss_case0@nareas <- 1
# Name of the case
ss_case0@Name <- "case0"
# Catch data
data_id <- as.numeric(names(sa_data$fishery$obs_total_catch_biomass$fleet1)) %in% model_year
ss_case0@Cat <- matrix(sa_data$fishery$obs_total_catch_biomass$fleet1[data_id],
                       nrow = 1
)
# State units of catch
ss_case0@Units <- sa_data$fishery$units_info$units_biomass
# Years
ss_case0@Year <- model_year
# Depletion relative to unfished
ss_case0@Dep <-  tail(ss_case0@Cat[1, ], n = 1) / max(ss_case0@Cat[1, ])

# VB maximum growth rate
ss_case0@vbK <- sa_data$biodata$k
# VB theoretical age at zero length
# default from EwE: -0.1
# https://www.researchgate.net/publication/267193103_Ecopath_with_Ecosim_A_User's_Guide#pf66
ss_case0@vbt0 <- sa_data$biodata$t0
# VB maximum length
# ss_case0@vbLinf <- (sa_data$biodata$winf * 1000 / sa_data$biodata$lw_a)^(1 / sa_data$biodata$lw_b)
ss_case0@vbLinf <- 500 # mm from FishBase
# Ratio of FMSY/M
# Q: Is it possible to find the ratio of FMSY to M
# for pelagic species in the Northwestern Atlantic Ocean?
# EwE FMSY values
### fmsy <- c(
###  0.092402786, 0.344614655, 1.101458788, 0.722086847,
###  1.569404483, 0.988227546, 0.72328496
###)

# Use compensation Fmsy (one value)
#fmsy_m <- mean(compensation_msy$FMSY / sa_data$biodata$natural_mortality_matrix[1, ])

# Use Fmsy at age
fmsy_m <- mean(om_fmsy_group / sa_data$biodata$natural_mortality_matrix[1, ])

ss_case0@FMSY_M <- fmsy_m

# ss_case@Mort <- mean(sa_data$biodata$natural_mortality_matrix[1,]) # = 1.1 and is greater than the max value of Mdb = 0.9 from DBSRA_()
ss_case0@Mort <- 0.9
ss_case0@CV_Mort <- 0.5 # default = 0.2

#ss_case0@Mort <- 0.67 #sa_data$biodata$natural_mortality_matrix[1, ]/(1-exp(-ss_case0@vbK*(1:7-ss_case0@vbt0)))^(-3*0.305)

# BMSY relative to unfished
# Dick and MacCall (2011): use 0.4 if target biomass is 40% of unfished biomass
### ss_case0@BMSY_B0 <- 0.3 # BAM FEC30% target

ss_case0@BMSY_B0 <- compensation_msy$BMSY / B0_F0_menhaden
# Length at 50% maturity
age_maturity50 <- 2
ss_case0@L50 <- ss_case0@vbLinf *
  (1 - exp(-ss_case0@vbK *
             (age_maturity50 - ss_case0@vbt0)))

# Run DBSRA, DBSRA_40 and DBSRA4010
## DBSRA: Base Version. TAC is calculated assumed MSY harvest rate
## multiplied by the estimated current abundance (estimated B0 x Depletion)

## DBSRA_40: Same as the Base Version but assumes 40 percent current
## depletion (Bcurrent/B0 = 0.4), which is more or less the most
## optimistic state for a stock (i.e. very close to BMSY/B0 for many stocks).

## DBSRA4010: Base version paired with the 40-10 rule that linearly
## throttles back the TAC when depletion is below 0.4 down to zero at
## 10 percent of unfished biomass.

ss_case01 <- ss_case0
dbsra <- DLMtool::DBSRA(1, ss_case01, plot = FALSE)
# dbsra40 <- DLMtool::DBSRA_40(1, ss_case0, plot = TRUE)
# dbsra4010 <- DLMtool::DBSRA4010(1, ss_case0, plot = TRUE)

# Projection: Scenario A --------------------------------------------------
projection_output1 <- list()
for (i in 1:length(projection_year)) {
  ss_case <- ss_case01
  year <- model_year[1]:(projection_year[i] - 1)
  data_id <- as.numeric(names(sa_data$fishery$obs_total_catch_biomass$fleet1)) %in% year
  ss_case@Cat <- matrix(sa_data$fishery$obs_total_catch_biomass$fleet1[data_id],
                        nrow = 1
  )
  ss_case@Year <- year
  #ss_case@Dep <- tail(ss_case@Cat[1, ], n = 1) / ss_case@Cat[1, 1]
  ss_case@Dep <- ss_case01@Dep
  # projection_output1[[i]] <- DLMtool::DBSRA_(
  #   x = 1,
  #   Data = ss_case,
  #   depo = NULL, # Optional fixed depletion (single value)
  #   hcr = NULL # Optional harvest control rule for throttling catch as a function of B/B0.
  # )
  projection_output1[[i]] <- DBSRA_return_FMSY(
    x = 1,
    Data = ss_case,
    depo = NULL, # Optional fixed depletion (single value)
    hcr = NULL # Optional harvest control rule for throttling catch as a function of B/B0.
  )
}
save(projection_output1, file = here::here("data", "data_poor", ewe_scenario_name, paste0("dbsra_projection_output1_", terminal_year, scenario_filename, ".RData")))

# plot figures
jpeg(filename = file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "dbsra_projection1.jpeg")), width = 250, height = 150, units = "mm", res = 1800)
par(mfrow=c(1,2))
ylim <- range(projection_output1[[1]]$C_hist, unlist(sapply(projection_output1, "[", "TAC")))
plot(c(model_year, projection_year),
     c(projection_output1[[1]]$C_hist, rep(NA, length(projection_year))),
     type = "l",
     ylim = ylim,
     xlab = "Year", ylab = "Catch (mt)"
)

for (i in 1:length(projection_year)) {
  boxplot(projection_output1[[i]]$TAC,
          add = TRUE,
          at = projection_year[i],
          col = "grey",
          width = 1,
          outline = TRUE,
          axes = FALSE
  )
}

catch_year <- projection_year[1:length(projection_year) - 1]
points(catch_year, sa_data$fishery$obs_total_catch_biomass$fleet1[names(sa_data$fishery$obs_total_catch_biomass$fleet1) %in% catch_year],
       pch = 16, col = "blue"
)

legend("top",
       c("Observed Catch (model year)", "Observed Catch (projection year)", "TAC"),
       lty = c(1, NA, NA),
       pch = c(NA, 16, NA),
       col = c("black", "blue", NA),
       fill = c(NA, NA, "gray"),
       border = c(NA, NA, "black"),
       bty = "n",
       title = paste0(ewe_scenario_name, "_OM: Scenario A")
)
legend("left",
       c(
         paste0("DBSRA@Dep = ", round(ss_case01@Dep, digits = 2)),
         paste0("DBSRA@BMSY_B0 = ", round(ss_case01@BMSY_B0, digits = 2)),
         paste0("DBSRA@FMSY_M = ", round(ss_case01@FMSY_M, digits = 2))
       ),
       bty="n")

# plot biomass
dbsraB_projection <- list()
for (i in 1:length(projection_year)){
  dbsraB_projection[[i]] <- apply(projection_output1[[i]]$Btrend, 2, median)
}

# ylim <- range(biomass_ewe[time_id], unlist(dbsraB_projection))
# for (i in 1:length(projection_year)){
#   lines(model_year[1]:(projection_year[i]-1), dbsraB_projection[[i]],
#         lty=i, col=i)
# }
# legend("bottomright",
#        c("EWE", paste0("DBSRA", projection_year)),
#        pch = c(16, rep(NA, length(projection_year))),
#        lty = c(NA, 1:length(projection_year)),
#        col = c(1, 1:length(projection_year)),
#        title = "Case 0: Scenario A", 
#        bty = "n"
# )
ylim <- range(biomass_ewe[time_id], unlist(dbsraB_projection[[1]]))
plot(c(model_year),
     biomass_ewe[time_id][1:length(model_year)],
     xlab = "Year", ylab = "Biomass (mt)",
     ylim = ylim,
     pch = 16
)
lines(model_year[1]:(projection_year[1]-1), dbsraB_projection[[1]],
        lty=1, col=1)
legend("bottomright",
       c("EWE", "DBSRA"),
       pch = c(16, NA),
       lty = c(NA, 1),
       col = c(1, 1),
       title = paste0(ewe_scenario_name, "_OM: Scenario A"), 
       bty = "n"
)
dev.off()

# Scenario B ---------------------------
ss_case02 <- ss_case0
ss_case02@Dep <- biomass_ewe[time_id][terminal_year-start_year+1]/B0_F0_menhaden

dbsra <- DLMtool::DBSRA(1, ss_case02, plot = FALSE)

# Projection: Scenario B --------------------------------------------------
projection_output2 <- list()
for (i in 1:length(projection_year)) {
  ss_case <- ss_case02
  year <- model_year[1]:(projection_year[i] - 1)
  data_id <- as.numeric(names(sa_data$fishery$obs_total_catch_biomass$fleet1)) %in% year
  ss_case@Cat <- matrix(sa_data$fishery$obs_total_catch_biomass$fleet1[data_id],
                        nrow = 1
  )
  ss_case@Year <- year
  #ss_case@Dep <- tail(ss_case@Cat[1, ], n = 1) / ss_case@Cat[1, 1]
  ss_case@Dep <- ss_case02@Dep
  # projection_output2[[i]] <- DLMtool::DBSRA_(
  #   x = 1,
  #   Data = ss_case,
  #   depo = NULL, # Optional fixed depletion (single value)
  #   hcr = NULL # Optional harvest control rule for throttling catch as a function of B/B0.
  # )
  projection_output2[[i]] <- DBSRA_return_FMSY(
    x = 1,
    Data = ss_case,
    depo = NULL, # Optional fixed depletion (single value)
    hcr = NULL # Optional harvest control rule for throttling catch as a function of B/B0.
  )
}

save(projection_output2, file = here::here("data", "data_poor", ewe_scenario_name, paste0("dbsra_projection_output2_", terminal_year, scenario_filename, ".RData")))

# plot figures
jpeg(filename = file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "dbsra_projection2.jpeg")), width = 250, height = 150, units = "mm", res = 1800)
par(mfrow=c(1,2))
ylim <- range(projection_output2[[1]]$C_hist, unlist(sapply(projection_output2, "[", "TAC")))
plot(c(model_year, projection_year),
     c(projection_output2[[1]]$C_hist, rep(NA, length(projection_year))),
     type = "l",
     ylim = ylim,
     xlab = "Year", ylab = "Catch (mt)"
)

for (i in 1:length(projection_year)) {
  boxplot(projection_output2[[i]]$TAC,
          add = TRUE,
          at = projection_year[i],
          col = "grey",
          width = 1,
          outline = TRUE,
          axes = FALSE
  )
}

catch_year <- projection_year[1:length(projection_year) - 1]
points(catch_year, sa_data$fishery$obs_total_catch_biomass$fleet1[names(sa_data$fishery$obs_total_catch_biomass$fleet1) %in% catch_year],
       pch = 16, col = "blue"
)

legend("top",
       c("Observed Catch (model year)", "Observed Catch (projection year)", "TAC"),
       lty = c(1, NA, NA),
       pch = c(NA, 16, NA),
       col = c("black", "blue", NA),
       fill = c(NA, NA, "gray"),
       border = c(NA, NA, "black"),
       bty = "n",
       title = paste0(ewe_scenario_name, "_OM: Scenario B")
)

legend("left",
       c(
         paste0("DBSRA@Dep = ", round(ss_case02@Dep, digits = 2)),
         paste0("DBSRA@BMSY_B0 = ", round(ss_case02@BMSY_B0, digits = 2)),
         paste0("DBSRA@FMSY_M = ", round(ss_case02@FMSY_M, digits = 2))
       ),
       bty="n")

# plot biomass
dbsraB_projection <- list()
for (i in 1:length(projection_year)){
  dbsraB_projection[[i]] <- apply(projection_output2[[i]]$Btrend, 2, median)
}

# ylim <- range(biomass_ewe[time_id], unlist(dbsraB_projection))
# plot(c(model_year, projection_year),
#      biomass_ewe[time_id],
#      xlab = "Year", ylab = "Biomass (mt)",
#      ylim = ylim,
#      pch = 16
# )
# for (i in 1:length(projection_year)){
#   lines(model_year[1]:(projection_year[i]-1), dbsraB_projection[[i]],
#         lty=i, col=i)
# }
# legend("bottomright",
#        c("EWE", paste0("DBSRA", projection_year)),
#        bty = "n",
#        pch = c(16, rep(NA, length(projection_year))),
#        lty = c(NA, 1:length(projection_year)),
#        col = c(1, 1:length(projection_year)),
#        title = "Case 0: Scenario B"
# )
ylim <- range(biomass_ewe[time_id], unlist(dbsraB_projection[[1]]))
plot(c(model_year),
     biomass_ewe[time_id][1:length(model_year)],
     xlab = "Year", ylab = "Biomass (mt)",
     ylim = ylim,
     pch = 16
)
lines(model_year[1]:(projection_year[1]-1), dbsraB_projection[[1]],
      lty=1, col=1)
legend("bottomright",
       c("EWE", "DBSRA"),
       pch = c(16, NA),
       lty = c(NA, 1),
       col = c(1, 1),
       title = paste0(ewe_scenario_name, "_OM: Scenario B"),
       bty = "n"
)
dev.off()

time_id_model_year <- seq(1, nrow(biomass[[1]]), by = 12)[1:length(model_year)]
biomass_ewe_model_year <- apply(biomass[[1]][, age_name], 1, sum)[time_id_model_year] * 1000000

ss_case03 <- ss_case0

equi_year <- start_year:(start_year+10)
catch_multiplier <- 1
equi_catch <- rep(sa_data$fishery$obs_total_catch_biomass$fleet1[1], length = length(equi_year)) * catch_multiplier
# equi_catch <- rep(max(sa_data$fishery$obs_total_catch_biomass$fleet), length = length(equi_year)) * catch_multiplier

data_id <- as.numeric(names(sa_data$fishery$obs_total_catch_biomass$fleet1)) %in% model_year
ss_case03@Cat <- matrix(c(equi_catch, sa_data$fishery$obs_total_catch_biomass$fleet1[data_id]),
                        nrow = 1
)

ss_case03@Year <- (model_year[1] - length(equi_year)):tail(model_year, n=1)

dep = seq(0.1, 0.9, by = 0.05)
bmsy_b0 <- seq(0.1, 0.9, by = 0.05)
fmsy_m <- seq(0.1, 2, by = 0.05)

param_matrix <- expand.grid(dep, bmsy_b0, fmsy_m)
colnames(param_matrix) <- c("dep", "bmsy_b0", "fmsy_m")
param_matrix$ssd <- NA

find_min_ssd = FALSE
print(terminal_year)
print(scenario_filename)
if (find_min_ssd == TRUE){
  for (i in 1:nrow(param_matrix)){
  ss_case <- ss_case03
  ss_case@Dep <- param_matrix$dep[i]
  ss_case@BMSY_B0 <- param_matrix$bmsy_b0[i]
  ss_case@FMSY_M <- param_matrix$fmsy_m[i]
  dbsra <- DLMtool::DBSRA_(
    x = 1,
    Data = ss_case,
    depo = NULL,
    hcr = NULL
  )
  dbsraB <- apply(dbsra$Btrend, 2, median)
  param_matrix$ssd[i] <- sum((dbsraB[(length(equi_year)+1):length(dbsraB)] - biomass_ewe_model_year)^2)
}

write.csv(param_matrix, here::here("data", "data_poor", ewe_scenario_name, paste0("dbsra_ssd_", terminal_year, scenario_filename, ".csv")))
}


param_matrix <- read.csv(here::here("data", "data_poor", ewe_scenario_name, paste0("dbsra_ssd_", terminal_year, scenario_filename, ".csv")))

ss_case <- ss_case03
ss_case@Dep <- param_matrix$dep[which.min(param_matrix$ssd)]
ss_case@BMSY_B0 <- param_matrix$bmsy_b0[which.min(param_matrix$ssd)]
ss_case@FMSY_M <- param_matrix$fmsy_m[which.min(param_matrix$ssd)] 

# dbsra <- DLMtool::DBSRA_(
#   x = 1,
#   Data = ss_case,
#   depo = NULL,
#   hcr = NULL
# )

dbsra <- DBSRA_return_FMSY(
  x = 1,
  Data = ss_case,
  depo = NULL,
  hcr = NULL
)

# dbsra_fmsym <- apply(dbsra$Btrend, 2, median) # Need to check

ss_case_final <- ss_case
dbsra <- DLMtool::DBSRA(1, ss_case_final, plot = FALSE)

# Projection: scenario C --------------------------------------------------
projection_output3 <- list()
equi_year_temp <- list()
for (i in 1:length(projection_year)) {
  ss_case_temp <- ss_case_final

  equi_year_temp[[i]] <- equi_year
  # equi_year_temp[[i]] <- 2000:(2012+i-1)
  equi_catch_temp <- rep(sa_data$fishery$obs_total_catch_biomass$fleet1[1], length = length(equi_year_temp[[i]])) * catch_multiplier
  # equi_catch_temp <- rep(max(sa_data$fishery$obs_total_catch_biomass$fleet1), length = length(equi_year_temp[[i]])) * catch_multiplier

  # equi_id_temp <- as.numeric(names(sa_data$fishery$obs_total_catch_biomass$fleet1)) %in% equi_year_temp[[i]]
  # equi_catch_temp <- sa_data$fishery$obs_total_catch_biomass$fleet1[equi_id_temp]

  data_id <- as.numeric(names(sa_data$fishery$obs_total_catch_biomass$fleet1)) %in% c(model_year[1]:(projection_year[i] - 1))

  ss_case_temp@Cat <- matrix(c(equi_catch_temp, sa_data$fishery$obs_total_catch_biomass$fleet1[data_id]),
                             nrow = 1
  )
  ss_case_temp@Year <- ss_case_temp@Year[1]:(projection_year[i]-1)


  # projection_output3[[i]] <- DLMtool::DBSRA_(
  #   x = 1,
  #   Data = ss_case_temp,
  #   depo = NULL, # Optional fixed depletion (single value)
  #   hcr = NULL # Optional harvest control rule for throttling catch as a function of B/B0.
  # )
  projection_output3[[i]] <- DBSRA_return_FMSY(
    x = 1,
    Data = ss_case_temp,
    depo = NULL, # Optional fixed depletion (single value)
    hcr = NULL # Optional harvest control rule for throttling catch as a function of B/B0.
  )
}

save(projection_output3, file = here::here("data", "data_poor", ewe_scenario_name, paste0("dbsra_projection_output3_", terminal_year, scenario_filename, ".RData")))

# plot figures
jpeg(filename = file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "dbsra_projection3.jpeg")), width = 250, height = 150, units = "mm", res = 1800)
par(mfrow=c(1,2))
c_hist <- projection_output3[[1]]$C_hist
ylim <- range(c_hist[(length(equi_year_temp[[1]])+1):length(c_hist)], unlist(sapply(projection_output3, "[", "TAC")))
plot(c(ss_case_temp@Year[(length(equi_year_temp[[1]])+1):length(c_hist)], projection_year),
     c(c_hist[(length(equi_year_temp[[1]])+1):length(c_hist)], rep(NA, length(projection_year))),
     type = "l",
     ylim = ylim,
     xlab = "Year", ylab = "Catch (mt)"
)

for (i in 1:length(projection_year)) {
  boxplot(projection_output3[[i]]$TAC,
          add = TRUE,
          at = projection_year[i],
          col = "grey",
          width = 1,
          outline = TRUE,
          axes = FALSE
  )
}

catch_year <- projection_year[1:length(projection_year) - 1]
points(catch_year, sa_data$fishery$obs_total_catch_biomass$fleet1[names(sa_data$fishery$obs_total_catch_biomass$fleet1) %in% catch_year],
       pch = 16, col = "blue"
)

legend("top",
       c("Observed Catch", "Observed Catch for TAC estimates", "TAC"),
       lty = c(1, NA, NA),
       pch = c(NA, 16, NA),
       col = c("black", "blue", NA),
       fill = c(NA, NA, "gray"),
       border = c(NA, NA, "black"),
       bty = "n",
       title = paste0(ewe_scenario_name, "_OM: Scenario C")
)

legend("left",
       c(
         paste0("DBSRA@Dep = ", round(param_matrix$dep[which.min(param_matrix$ssd)], digits = 2)),
         paste0("DBSRA@BMSY_B0 = ", round(param_matrix$bmsy_b0[which.min(param_matrix$ssd)], digits = 2)),
         paste0("DBSRA@FMSY_M = ", round(param_matrix$fmsy_m[which.min(param_matrix$ssd)], digits = 2))),
       bty="n")

# plot biomass
dbsraB_projection <- list()
for (i in 1:length(projection_year)){
  full_dbsraB <- projection_output3[[i]]$Btrend
  year_id <- projection_output3[[i]]
  dbsraB_projection[[i]] <- apply(full_dbsraB[, (length(equi_year_temp[[i]])+1):ncol(full_dbsraB)], 2, median)
}

# ylim = range(biomass_ewe[time_id], unlist(dbsraB_projection))
# plot(c(model_year, projection_year),
#      biomass_ewe[time_id],
#      xlab = "Year", ylab = "Biomass (mt)",
#      ylim = ylim,
#      pch = 16
# )
# for (i in 1:length(projection_year)){
#   lines(model_year[1]:(projection_year[i]-1), dbsraB_projection[[i]],
#         lty=i, col=i)
# }
# legend("bottomright",
#        c("EWE", paste0("DBSRA", projection_year)),
#        bty = "n",
#        pch = c(16, rep(NA, length(projection_year))),
#        lty = c(NA, 1:length(projection_year)),
#        col = c(1, 1:length(projection_year)),
#        title = "Case 0: scenario C"
# )
ylim <- range(biomass_ewe[time_id], unlist(dbsraB_projection[[1]]))
plot(c(model_year),
     biomass_ewe[time_id][1:length(model_year)],
     xlab = "Year", ylab = "Biomass (mt)",
     ylim = ylim,
     pch = 16
)
lines(model_year[1]:(projection_year[1]-1), dbsraB_projection[[1]],
      lty=1, col=1)
legend("bottomright",
       c("EWE", "DBSRA"),
       bty = "n",
       pch = c(16, NA),
       lty = c(NA, 1),
       col = c(1, 1),
       title = paste0(ewe_scenario_name, "_OM: Scenario C")
)
dev.off()

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

Status of indicators (SOI)

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

where $i$ represents individual iterations

Cases 1-5 are based on the settings from Scenario A

# Projection based on Scenario A ------------------------------------

output <- compare_projections_dbsra(
  case_name = "scenario A",
  projection_output_data = projection_output1,
  projection_year = projection_year,
  model_year = model_year,
  lag = 1,
  amo = amo_unsmooth_lag1,
  pdsi = palmer_drought_severity_index,
  bassB = bass_bio,
  menhadenC = menhadenC, 
  menhadenE = menhadenE, 
  menhadenCPUE = menhadenCPUE,
  bassCPUE = bassCPUE,
  herringCPUE = herringCPUE, 
  menhadenV = menhadenV, 
  indices = c("amo", "pdsi", "bassB", "menhadenC", "menhadenE", "menhadenCPUE", "bassCPUE", "herringCPUE", "menhadenV")
)

soi_output1 <- output
soi_output1$lm_data_om <- lm_data_om
save(soi_output1, file = here::here("data", "data_poor", ewe_scenario_name, paste0("dbsra_soi_output1_", terminal_year, scenario_filename, ".RData")))

# Linear regression figure
lm_data <- rbind(lm_data_om, output$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"))

subset_data <- lm_data[!(lm_data$Variable %in% c("Menhaden fishing effort", "Menhaden CPUE")), ]
ggplot(subset_data, aes(x = Index_Value, y = Menhaden_Biomass, color = model)) +
  geom_point() +
  geom_smooth(method = lm) +
  facet_wrap(~Variable, scales = "free_x") +
  labs(
    # title = "Linear regression analysis",
    x = "Indicator Value",
    y = "Menhaden-like species biomass (mt)"
  ) +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=1)) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "linear_regression1.jpeg")))

soi_data_melt <- output$soi_data_melt
subset_data <- soi_data_melt[!(soi_data_melt$variable %in% c("menhadenE", "menhadenCPUE")), ]
ggplot(subset_data[subset_data$projection_year_id==projection_year[1], ], aes(x = year, y = value)) +
  geom_line() +
  facet_wrap(~variable, scales = "free_y") +
  labs(
    title = "Scenario A: Status of Indicators",
    x = "Year",
    y = "Status of Indicators"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "soi1.jpeg")))

Bt_BMSY_dataframe <- as.data.frame(do.call(cbind, output$Bt_BMSY_list))
colnames(Bt_BMSY_dataframe) <- projection_year - 1
Bt_BMSY_stack <- stack(Bt_BMSY_dataframe)
ggplot(Bt_BMSY_stack, aes(x = ind, y = values)) +
  geom_boxplot() +
  labs(
    # title = "Scenario A: Bt/BMSY",
    x = "",
    y = "Bratio"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "bratio1.jpeg")))

tac_data_melt <- output$tac_data_melt
subset_data <- tac_data_melt[!(tac_data_melt$variable %in% c("menhadenE", "menhadenCPUE")), ]
ggplot(subset_data, aes(x=variable, y = value, color = projection_year_id)) +
  geom_boxplot() +
  labs(
    title = "Scenario A: TAC",
    x = "",
    y = "TAC"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "tac1_boxplot.jpeg")))

ggplot(subset_data, aes(x=variable, y = value, color = projection_year_id)) +
  geom_violin()+
  geom_boxplot(width=0.1) +
  # facet_wrap(~variable)+
  labs(
    title = "Scenario A: TAC",
    x = "",
    y = "TAC"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "tac1_violin.jpeg")))

tac_data_melt_median <- output$tac_data_melt_median
subset_data <- tac_data_melt_median[!(tac_data_melt_median$variable %in% c("menhadenE", "menhadenCPUE")), ]
ggplot(subset_data, aes(x = projection_year_id, y = value, color = variable)) +
  geom_line() +
  labs(
    title = "Scenario A: Median TAC from DBSRA and adjusted TAC",
    x = "",
    y = "TAC"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "tac1_line_median.jpeg")))

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

Cases 1-5 are based on the settings from Scenario B

# Projection based on Scenario B ------------------------------------

output <- compare_projections_dbsra(
  case_name = "scenario B",
  projection_output_data = projection_output2,
  projection_year = projection_year,
  model_year = model_year,
  lag = 1,
  amo = amo_unsmooth_lag1,
  pdsi = palmer_drought_severity_index,
  bassB = bass_bio,
  menhadenC = menhadenC, 
  menhadenE = menhadenE, 
  menhadenCPUE = menhadenCPUE,
  bassCPUE = bassCPUE,
  herringCPUE = herringCPUE, 
  menhadenV = menhadenV, 
  indices = c("amo", "pdsi", "bassB", "menhadenC", "menhadenE", "menhadenCPUE", "bassCPUE", "herringCPUE", "menhadenV")
)

soi_output2 <- output
soi_output2$lm_data_om <- lm_data_om
save(soi_output2, file = here::here("data", "data_poor", ewe_scenario_name, paste0("dbsra_soi_output2_", terminal_year, scenario_filename, ".RData")))

# Linear regression figure
lm_data <- rbind(lm_data_om, output$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"))

subset_data <- lm_data[!(lm_data$Variable %in% c("Menhaden fishing effort", "Menhaden CPUE")), ]
ggplot(subset_data, aes(x = Index_Value, y = Menhaden_Biomass, color = model)) +
  geom_point() +
  geom_smooth(method = lm) +
  facet_wrap(~Variable, scales = "free_x") +
  labs(
    # title = "Linear regression analysis",
    x = "Indicator Value",
    y = "Menhaden-like species biomass (mt)"
  ) +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=1)) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "linear_regression2.jpeg")))

soi_data_melt <- output$soi_data_melt
subset_data <- soi_data_melt[!(soi_data_melt$variable %in% c("menhadenE", "menhadenCPUE")), ]
ggplot(subset_data[subset_data$projection_year_id==projection_year[1], ], aes(x = year, y = value)) +
  geom_line() +
  facet_wrap(~variable, scales = "free_y") +
  labs(
    title = "Scenario B: Status of Indicators",
    x = "Year",
    y = "Status of Indicators"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "soi2.jpeg")))

Bt_BMSY_dataframe <- as.data.frame(do.call(cbind, output$Bt_BMSY_list))
colnames(Bt_BMSY_dataframe) <- projection_year - 1
Bt_BMSY_stack <- stack(Bt_BMSY_dataframe)
ggplot(Bt_BMSY_stack, aes(x = ind, y = values)) +
  geom_boxplot() +
  labs(
    # title = "Scenario B: Bt/BMSY",
    x = "",
    y = "Bratio"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "bratio2.jpeg")))

tac_data_melt <- output$tac_data_melt
subset_data <- tac_data_melt[!(tac_data_melt$variable %in% c("menhadenE", "menhadenCPUE")), ]
ggplot(subset_data, aes(x=variable, y = value, color = projection_year_id)) +
  geom_boxplot() +
  labs(
    title = "Scenario B: TAC",
    x = "",
    y = "TAC"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "tac2_boxplot.jpeg")))

ggplot(subset_data, aes(x=variable, y = value, color = projection_year_id)) +
  geom_violin()+
  geom_boxplot(width=0.1) +
  # facet_wrap(~variable)+
  labs(
    title = "Scenario B: TAC",
    x = "",
    y = "TAC"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "tac2_violin.jpeg")))

tac_data_melt_median <- output$tac_data_melt_median
subset_data <- tac_data_melt_median[!(tac_data_melt_median$variable %in% c("menhadenE", "menhadenCPUE")), ]
ggplot(subset_data, aes(x = projection_year_id, y = value, color = variable)) +
  geom_line() +
  labs(
    title = "Scenario B: Median TAC from DBSRA and adjusted TAC",
    x = "",
    y = "TAC"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "tac2_line_median.jpeg")))

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

Cases 1-5 are based on the settings from Scenario C

# Projection based on Scenario C ------------------------------------
projection_output3_clean <- list()

for (i in 1:length(projection_year)){
  temp <- projection_output3[[i]]
  temp$Btrend <- temp$Btrend[, (length(equi_year_temp[[i]])+1): ncol(temp$Btrend)]
  projection_output3_clean[[i]] <- temp
}

output <- compare_projections_dbsra(
  case_name = "scenario C",
  projection_output_data = projection_output3_clean,
  projection_year = projection_year,
  model_year = model_year,
  lag = 1,
  amo = amo_unsmooth_lag1,
  pdsi = palmer_drought_severity_index,
  bassB = bass_bio,
  menhadenC = menhadenC, 
  menhadenE = menhadenE, 
  menhadenCPUE = menhadenCPUE,
  bassCPUE = bassCPUE,
  herringCPUE = herringCPUE, 
  menhadenV = menhadenV, 
  indices = c("amo", "pdsi", "bassB", "menhadenC", "menhadenE", "menhadenCPUE", "bassCPUE", "herringCPUE", "menhadenV")
)

soi_output3 <- output
soi_output3$lm_data_om <- lm_data_om
save(soi_output3, file = here::here("data", "data_poor", ewe_scenario_name, paste0("dbsra_soi_output3_", terminal_year, scenario_filename, ".RData")))

# Linear regression figure
lm_data <- rbind(lm_data_om, output$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"))

subset_data <- lm_data[!(lm_data$Variable %in% c("Menhaden fishing effort", "Menhaden CPUE")), ]
ggplot(subset_data, aes(x = Index_Value, y = Menhaden_Biomass, color = model)) +
  geom_point() +
  geom_smooth(method = lm) +
  facet_wrap(~Variable, scales = "free_x") +
  labs(
    # title = "Linear regression analysis",
    x = "Indicator Value",
    y = "Menhaden-like species biomass (mt)"
  ) +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=1)) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "linear_regression3.jpeg")))

soi_data_melt <- output$soi_data_melt
subset_data <- soi_data_melt[!(soi_data_melt$variable %in% c("menhadenE", "menhadenCPUE")), ]
ggplot(subset_data[subset_data$projection_year_id==projection_year[1], ], aes(x = year, y = value)) +
  geom_line() +
  facet_wrap(~variable, scales = "free_y") +
  labs(
    title = "Scenario C: Status of Indicators",
    x = "Year",
    y = "Status of Indicators"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "soi3.jpeg")))

Bt_BMSY_dataframe <- as.data.frame(do.call(cbind, output$Bt_BMSY_list))
colnames(Bt_BMSY_dataframe) <- projection_year - 1
Bt_BMSY_stack <- stack(Bt_BMSY_dataframe)
ggplot(Bt_BMSY_stack, aes(x = ind, y = values)) +
  geom_boxplot() +
  labs(
    # title = "Scenario C: Bt/BMSY",
    x = "",
    y = "Bratio"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "bratio3.jpeg")))

tac_data_melt <- output$tac_data_melt
subset_data <- tac_data_melt[!(tac_data_melt$variable %in% c("menhadenE", "menhadenCPUE")), ]
ggplot(subset_data, aes(x=variable, y = value, color = projection_year_id)) +
  geom_boxplot() +
  labs(
    title = "Scenario C: TAC",
    x = "",
    y = "TAC"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "tac3_boxplot.jpeg")))

ggplot(subset_data, aes(x=variable, y = value, color = projection_year_id)) +
  geom_violin()+
  geom_boxplot(width=0.1) +
  # facet_wrap(~variable)+
  labs(
    title = "Scenario C: TAC",
    x = "",
    y = "TAC"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "tac3_violin.jpeg")))

tac_data_melt_median <- output$tac_data_melt_median
subset_data <- tac_data_melt_median[!(tac_data_melt_median$variable %in% c("menhadenE", "menhadenCPUE")), ]
ggplot(subset_data, aes(x = projection_year_id, y = value, color = variable)) +
  geom_line() +
  labs(
    title = "Scenario C: Median TAC from DBSRA and adjusted TAC",
    x = "",
    y = "TAC"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "tac3_line_median.jpeg")))

for (i in 1:nrow(tac_data_melt)){
  temp <- projection_output3[[which(projection_year == tac_data_melt$projection_year_id[i])]]
  temp_B <- (temp$Btrend[, 1] * temp$Bt_Kstore)[tac_data_melt$iter[i]]
  # tac_data_melt$fmort[i] <- tac_data_melt$value[i]/temp_B
  tac_data_melt$fmort[i] <- output$fmsy_data_melt$value[i]
}

fmort_data_melt_median <- aggregate(fmort ~ projection_year_id+variable, data = tac_data_melt, median)
# projection_f <- fmort_data_melt_median

projection_f <- rbind(
  fmort_data_melt_median,
  data.frame(
    projection_year_id = projection_year,
    variable = "om",
    fmort = 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$fmort[i] * fishery_sel
}

write.csv(projection_f, here::here("data", "data_poor", ewe_scenario_name, paste0("fmsy_median_data3_", terminal_year, scenario_filename, ".csv")))

Bt_BMSY_matrix <- matrix(NA, ncol=length(projection_year), nrow=length(projection_output3_clean[[1]]$Bt_Kstore))
for (projection_year_id in 1:length(projection_year)){
    projection_output <- projection_output3_clean
    Bt_BMSY_matrix[, projection_year_id] <- projection_output[[projection_year_id]]$Bt_Kstore / projection_output[[projection_year_id]]$BMSY_K_Mstore
}
write.csv(Bt_BMSY_matrix, here::here("data", "data_poor", ewe_scenario_name, paste0("Bt_BMSY_data3_", terminal_year, scenario_filename, ".csv")))

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

# tac_data_melt$biomass <- tac_data_melt$catch <- tac_data_melt$Bendyear_K <-
  # tac_data_melt$Bendyear_BMSY <- tac_data_melt$Bonanza_percentage <- tac_data_melt$Collapse_percentage <- NA

cases <- unique(tac_data_melt$variable)
iters <- unique(tac_data_melt$iter)

for (case_id in 1:length(cases)){
  temp <- tac_data_melt[tac_data_melt$variable==cases[case_id], ]

  for (iter in 1:length(iters)){
    tac_temp <- temp[temp$iter==iters[iter], ]

    ss_temp <- ss_case_final
    ss_temp@Cat <- matrix(c(ss_case_final@Cat, tac_temp$value), nrow = 1)
    ss_temp@Year <- c(ss_case_final@Year, as.numeric(as.character(tac_temp$projection_year_id)))
    model <- DLMtool::DBSRA_(
      x = 1,
      Data = ss_temp,
      depo = NULL, # Optional fixed depletion (single value)
      hcr = NULL # Optional harvest control rule for throttling catch as a function of B/B0.
    )

    tac_data_melt$biomass[(tac_data_melt$variable==cases[case_id]) & (tac_data_melt$iter == iters[iter])] <-
      apply(model$Btrend[, ss_temp@Year %in% as.numeric(as.character(tac_temp$projection_year_id))], 2, median)

    tac_data_melt$catch[(tac_data_melt$variable==cases[case_id]) & (tac_data_melt$iter == iters[iter])] <-
      tac_temp$value

    tac_data_melt$Bendyear_K[(tac_data_melt$variable==cases[case_id]) & (tac_data_melt$iter == iters[iter])] <-
      median(model$Bt_Kstore)

    tac_data_melt$Bendyear_BMSY[(tac_data_melt$variable==cases[case_id]) & (tac_data_melt$iter == iters[iter])] <-
      median(model$Bt_Kstore/model$BMSY_K_Mstore)

    tac_data_melt$K[(tac_data_melt$variable==cases[case_id]) & (tac_data_melt$iter == iters[iter])] <-
      median(model$Btrend[, ncol(model$Btrend)] / model$Bt_Kstore)

    tac_data_melt$Bt_K[(tac_data_melt$variable==cases[case_id]) & (tac_data_melt$iter == iters[iter])] <-
      tac_data_melt$biomass[(tac_data_melt$variable==cases[case_id]) & (tac_data_melt$iter == iters[iter])] /
      tac_data_melt$K[(tac_data_melt$variable==cases[case_id]) & (tac_data_melt$iter == iters[iter])]


    # tac_data_melt$Bonanza_Percentage[(tac_data_melt$variable==cases[case_id]) & (tac_data_melt$iter == iters[iter])] <-
    #   length(model$Bt_Kstore[model$Bt_Kstore>0.8])/length(model$Bt_Kstore)*100
    #
    # tac_data_melt$Collapse_Percentage[(tac_data_melt$variable==cases[case_id]) & (tac_data_melt$iter == iters[iter])] <-
    #   length(model$Bt_Kstore[model$Bt_Kstore<0.2])/length(model$Bt_Kstore)*100


  }



}

ggplot(tac_data_melt, aes(x=projection_year_id, y = biomass, color = projection_year_id)) +
  geom_violin()+
  geom_boxplot(width=0.1) +
  facet_wrap(~variable)+
  labs(
    # title = " Estimated biomass based on TAC",
    x = "",
    y = "Biomass"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "assessment_biomass_projections_fixy.jpeg")))

ggplot(tac_data_melt, aes(x=projection_year_id, y = biomass, color = projection_year_id)) +
  geom_violin()+
  geom_boxplot(width=0.1) +
  facet_wrap(~variable, scales="free_y")+
  labs(
    # title = " Estimated biomass based on TAC",
    x = "",
    y = "Biomass"
  ) +
  theme_bw()
ggsave(file.path(figure_path, paste0(ewe_scenario_name, "_", terminal_year, scenario_filename, "assessment_biomass_projections_freey.jpeg")))

om_K <- B0_F0_menhaden
om_bmsy <- compensation_msy$BMSY
om_fmsy <- compensation_msy$FMSY

cases <- c("om", "dbsra",
           "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_poor_", 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_poor_", 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(tac_data_melt$variable)
for (i in 1:length(case_names)) {
  temp <- tac_data_melt[tac_data_melt$variable == case_names[i], ]

  evaluation_table <- rbind(
    evaluation_table,
    c(
      mean(aggregate(catch ~ projection_year_id, data = temp, median)$catch),
      aggregate(Bendyear_K ~ projection_year_id, data = temp[temp$projection_year_id == tail(projection_year, n = 1), ], median)$Bendyear_K,
      aggregate(Bendyear_BMSY ~ projection_year_id, data = temp[temp$projection_year_id == tail(projection_year, n = 1), ], median)$Bendyear_BMSY,
      mean(aggregate(biomass ~ projection_year_id, data = temp, median)$biomass),
      sum(aggregate(Bt_K ~ projection_year_id, data = temp, median)$Bt_K > 0.8),
      sum(aggregate(Bt_K ~ projection_year_id, data = temp, median)$Bt_K < 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_poor", 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.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close