In Bai-Li-NOAA/IFA4EBFM: Interdisciplinary Projections for Ecosystem-Based Fisheries Management
knitr::opts_chunk$set(
echo = FALSE,
message = FALSE,
warning = FALSE,
fig.width = 12,
fig.height = 10,
results = "asis"
)
options(scipen = 1, digits = 4)
devtools::load_all()
library(tidyr)
Evaluation of equilibrium MSY reference points from EwE user manual
- Single-species FMSY (stationary FMSY)
- Run the Ecosim model to equilibrium for a range of F values while holding biomasses of all other groups constant (except other life history stanzas for the same species in split and multi-stanza cases)
- Examine predicted compensatory responses by the group caused by the foraging arena functional response and related foraging time adjustment parameters
- Ecological FMSY (compensatory FMSY)
- Allows all groups' biomass to change in response to the biomass changes of the focal FMSY groups
- The assessment usually produces a more optimistic estimate of MSY for each group
- Increases in food species abundances and decreases in predator abundances when any one group is harvested harder
Output comparison
Run model to get compensation FMSY vs stationary FMSY
# load sa_data to get maturity data
scenario_filename <- "_lowCV"
terminal_year <- 2012
fishery_CV_input <- 0.05
survey_CV_input <- 0.1
add_environmental_effects <- FALSE
add_fleet_dynamics <- FALSE
source(here::here("Rscript", "simulationOM3.R"))
maturity_matrix <- sa_data$biodata$maturity_matrix
# load 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)
# base scenario
compensation_fmsy_s1 <- read_fmsy_table(
file_names = here::here(
"data", "ewe", "7ages_ecosim_om_msy",
paste0("msy_base_run_age", ages), "FMSY_fullCompensation.csv"
),
colnames = "Fmsy",
functional_groups = functional_groups
)
compensation_fmsy_s1 <- c(
compensation_fmsy_s1,
read_fmsy_table(
file_names = here::here(
"data", "ewe", "7ages_ecosim_om_msy",
"msy_base_run_fleet", "FMSY_fullCompensation.csv"
),
colnames = "Fmsy",
functional_groups = functional_groups
)
)
stationary_fmsy_s1 <- read_fmsy_table(
file_names = here::here(
"data", "ewe", "7ages_ecosim_om_msy",
paste0("msy_base_run_age", ages), "FMSY_StationarySystem.csv"
),
colnames = "Fmsy",
functional_groups = functional_groups
)
stationary_fmsy_s1 <- c(
stationary_fmsy_s1,
read_fmsy_table(
file_names = here::here(
"data", "ewe", "7ages_ecosim_om_msy",
"msy_base_run_fleet", "FMSY_StationarySystem.csv"
),
colnames = "Fmsy",
functional_groups = functional_groups
)
)
mean_difference <- all.equal(current = stationary_fmsy_s1, target = compensation_fmsy_s1)
mean_difference_value <- as.numeric(sapply(strsplit(mean_difference, " "), "[[", 6))
var_s1 <- var(mean_difference_value)
# Environmental forcing scenario
compensation_fmsy_s2 <- read_fmsy_table(
file_names = here::here(
"data", "ewe", "7ages_ecosim_om_msy",
paste0("msy_forcing_pdsi_egg_amo1_age", ages), "FMSY_fullCompensation.csv"
),
colnames = "Fmsy",
functional_groups = functional_groups
)
compensation_fmsy_s2 <- c(
compensation_fmsy_s2,
read_fmsy_table(
file_names = here::here(
"data", "ewe", "7ages_ecosim_om_msy",
"msy_forcing_pdsi_egg_amo1_fleet", "FMSY_fullCompensation.csv"
),
colnames = "Fmsy",
functional_groups = functional_groups
)
)
stationary_fmsy_s2 <- read_fmsy_table(
file_names = here::here(
"data", "ewe", "7ages_ecosim_om_msy",
paste0("msy_forcing_pdsi_egg_amo1_age", ages), "FMSY_StationarySystem.csv"
),
colnames = "Fmsy",
functional_groups = functional_groups
)
stationary_fmsy_s2 <- c(
stationary_fmsy_s2,
read_fmsy_table(
file_names = here::here(
"data", "ewe", "7ages_ecosim_om_msy",
"msy_forcing_pdsi_egg_amo1_fleet", "FMSY_StationarySystem.csv"
),
colnames = "Fmsy",
functional_groups = functional_groups
)
)
mean_difference <- all.equal(current = stationary_fmsy_s2, target = compensation_fmsy_s2)
mean_difference_value <- as.numeric(sapply(strsplit(mean_difference, " "), "[[", 6))
var_s2 <- var(mean_difference_value)
if (var_s1 == 0 & var_s2 == 0) {
cat("- Mean relative differences between compensation and stationary FMSYs across 7 age classes from scenarios S1 and S2 = ", unique(mean_difference_value), "\n")
} else {
cat("- Mean relative differences between compensation and stationary FMSYs vary by age groups\n")
}
Run model by functional group vs by fleet
mean_difference <- c()
list_length <- length(compensation_fmsy_s1)
# base scenario
for (i in seq_along(compensation_fmsy_s1)) {
mean_difference[i] <- all.equal(current = compensation_fmsy_s1[[list_length]], target = compensation_fmsy_s1[[i]])
}
compensation_mean_difference_s1 <- all(mean_difference)
for (i in seq_along(stationary_fmsy_s1)) {
mean_difference[i] <- all.equal(current = stationary_fmsy_s1[[list_length]], target = stationary_fmsy_s1[[i]])
}
stationary_mean_difference_s1 <- all(mean_difference)
# environmental forcing scenario
for (i in seq_along(compensation_fmsy_s2)) {
mean_difference[i] <- all.equal(current = compensation_fmsy_s2[[list_length]], target = compensation_fmsy_s2[[i]])
}
compensation_mean_difference_s2 <- all(mean_difference)
for (i in seq_along(stationary_fmsy_s2)) {
mean_difference[i] <- all.equal(current = stationary_fmsy_s2[[list_length]], target = stationary_fmsy_s2[[i]])
}
stationary_mean_difference_s2 <- all(mean_difference)
if (all(c(compensation_mean_difference_s1, stationary_mean_difference_s1, compensation_mean_difference_s2, stationary_mean_difference_s2))) {
cat("- FMSYs by functional group are identical to FMSYs by fleet\n")
} else {
cat("- FMSYs by functional group are different to FMSYs by fleet\n")
}
Run model to equilibrium for a range of F values vs run to depletion
msy_ages_label <- c(0:5, "6+")
maxF_file_path <- here::here(
"data", "ewe", "7ages_ecosim_om_msy",
paste0("msy_forcing_pdsi_egg_amo1_age", ages, "_details")
)
maxF_file_path <- c(
maxF_file_path,
here::here("data", "ewe", "7ages_ecosim_om_msy", "msy_forcing_pdsi_egg_amo1_fleet_details")
)
catch_file_names <- c(paste0("menhaden ", msy_ages_label, "_Catch_FullCompensation.csv"), "F.menhaden_Catch_FullCompensation.csv")
biomass_file_names <- c(paste0("menhaden ", msy_ages_label, "_B_FullCompensation.csv"), "F.menhaden_B_FullCompensation.csv")
compensation_msy_maxF <- list()
for (i in seq_along(maxF_file_path)) {
compensation_msy_maxF[[i]] <- read_ewe_reference_points(
file_path = maxF_file_path[i],
file_names = c(
catch_file_names[i],
biomass_file_names[i]
),
reference_points_scenario = "compensation",
functional_groups = functional_groups,
key_functional_group = "AtlanticMenhaden",
ages = ages,
plot = FALSE,
maturity_vec = maturity_matrix[1,]
)
}
depletion_file_path <- here::here(
"data", "ewe", "7ages_ecosim_om_msy",
paste0("msy_forcing_pdsi_egg_amo1_age", ages, "_details_depletion")
)
depletion_file_path <- c(
depletion_file_path,
here::here("data", "ewe", "7ages_ecosim_om_msy", "msy_forcing_pdsi_egg_amo1_fleet_details_depletion")
)
compensation_msy_depletion <- list()
for (i in seq_along(maxF_file_path)) {
compensation_msy_depletion[[i]] <- read_ewe_reference_points(
file_path = depletion_file_path[i],
file_names = c(
catch_file_names[i],
biomass_file_names[i]
),
reference_points_scenario = "compensation",
functional_groups = functional_groups,
key_functional_group = "AtlanticMenhaden",
ages = ages,
plot = FALSE,
maturity_vec = maturity_matrix[1,]
)
}
result_compensation <- all.equal(current = compensation_msy_maxF, target = compensation_msy_depletion)
catch_file_names <- c(paste0("menhaden ", msy_ages_label, "_Catch_StationarySystem.csv"), "F.menhaden_Catch_StationarySystem.csv")
biomass_file_names <- c(paste0("menhaden ", msy_ages_label, "_B_StationarySystem.csv"), "F.menhaden_B_StationarySystem.csv")
stationary_msy_maxF <- list()
for (i in seq_along(maxF_file_path)) {
stationary_msy_maxF[[i]] <- read_ewe_reference_points(
file_path = maxF_file_path[i],
file_names = c(
catch_file_names[i],
biomass_file_names[i]
),
reference_points_scenario = "stationary",
functional_groups = functional_groups,
key_functional_group = "AtlanticMenhaden",
ages = ages,
plot = FALSE,
maturity_vec = maturity_matrix[1,]
)
}
stationary_msy_depletion <- list()
for (i in seq_along(maxF_file_path)) {
stationary_msy_depletion[[i]] <- read_ewe_reference_points(
file_path = depletion_file_path[i],
file_names = c(
catch_file_names[i],
biomass_file_names[i]
),
reference_points_scenario = "stationary",
functional_groups = functional_groups,
key_functional_group = "AtlanticMenhaden",
ages = ages,
plot = FALSE,
maturity_vec = maturity_matrix[1,]
)
}
result_stationary <- all.equal(current = stationary_msy_maxF, target = stationary_msy_depletion)
if (all(c(result_compensation, result_stationary))) {
cat("- FMSYs from the two searching methods are identical\n")
} else {
cat("- FMSYs from the two searching methods are different\n")
}
Single FMSY vs multiple FMSYs
- EwE default summary output files (e.g., FMSY_FullCompensation.csv or FMSY_StationarySystem.csv) show the estimate of FMSY for each functional group (Table 1-4).
- EwE also produces the catch of each functional group over a range of F values in F.menhaden_Catch_FullCompensation.csv or F.menhaden_Catch_StationarySystem.csv. To extract reference points, we can get a vector of fishing effort multipliers that correspond to the maximum catch of the case stock, and them multiply the fishing effort multipliers with seven base fishing mortality rates from 1985 to obtain seven FMSY values by age.
compensation_fmsy_s1 <- read_fmsy_table(
file_names = here::here(
"data", "ewe", "7ages_ecosim_om_msy",
"msy_base_run_fleet", "FMSY_fullCompensation.csv"
),
colnames = c(
"Group", "TL", "Fbase", "Cbase",
"Vbase", "FmsyFound", "Fmsy", "Cmsy", "Vmsy"
),
functional_groups = functional_groups
)
stationary_fmsy_s1 <- read_fmsy_table(
file_names = here::here(
"data", "ewe", "7ages_ecosim_om_msy",
"msy_base_run_fleet", "FMSY_StationarySystem.csv"
),
colnames = c(
"Group", "TL", "Fbase", "Cbase",
"Vbase", "FmsyFound", "Fmsy", "Cmsy", "Vmsy"
),
functional_groups = functional_groups
)
knitr::kable(compensation_fmsy_s1,
caption = "Table 1. Base scenario S1 compensation FMSY output from EwE"
)
knitr::kable(stationary_fmsy_s1,
caption = "Table 2. Base scenario S1 stationary FMSY output from EwE"
)
compensation_fmsy_s2 <- read_fmsy_table(
file_names = here::here(
"data", "ewe", "7ages_ecosim_om_msy",
"msy_forcing_pdsi_egg_amo1_fleet", "FMSY_fullCompensation.csv"
),
colnames = c(
"Group", "TL", "Fbase", "Cbase",
"Vbase", "FmsyFound", "Fmsy", "Cmsy", "Vmsy"
),
functional_groups = functional_groups
)
stationary_fmsy_s2 <- read_fmsy_table(
file_names = here::here(
"data", "ewe", "7ages_ecosim_om_msy",
"msy_forcing_pdsi_egg_amo1_fleet", "FMSY_StationarySystem.csv"
),
colnames = c(
"Group", "TL", "Fbase", "Cbase",
"Vbase", "FmsyFound", "Fmsy", "Cmsy", "Vmsy"
),
functional_groups = functional_groups
)
knitr::kable(compensation_fmsy_s2,
caption = "Table 3. Environmental forcing scenario S2 compensation FMSY output from EwE"
)
knitr::kable(stationary_fmsy_s2,
caption = "Table 4. Environmental forcing scenario S2 stationary FMSY output from EwE"
)
msy_file_path <- here::here("data", "ewe", "7ages_ecosim_om_msy", "msy_base_run_fleet_details")
compensation_msy_s1_MaxSumCatch <- 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,
maturity_vec = maturity_matrix[1,]
)
compensation_fmsy_s1_MaxSumCatch <- compensation_msy_s1_MaxSumCatch$FMSY * compensation_fmsy_s1[[1]]$Fbase
compensation_msy_s1_MaxCatchVec <- read_ewe_reference_points_vec(
file_path = msy_file_path,
file_names = c(
"F.menhaden_Catch_FullCompensation.csv",
"F.menhaden_B_FullCompensation.csv"
),
reference_points_scenario = "compensation",
functional_groups = functional_groups,
key_functional_group = "AtlanticMenhaden",
ages = ages,
maturity_vec = maturity_matrix[1,]
)
compensation_fmsy_s1_CatchVec <- compensation_msy_s1_MaxCatchVec$FMSY * compensation_fmsy_s1[[1]]$Fbase
stationary_msy_s1_MaxSumCatch <- 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,
maturity_vec = maturity_matrix[1,]
)
stationary_fmsy_s1_MaxSumCatch <- stationary_msy_s1_MaxSumCatch$FMSY * stationary_fmsy_s1[[1]]$Fbase
stationary_msy_s1_MaxCatchVec <- read_ewe_reference_points_vec(
file_path = msy_file_path,
file_names = c(
"F.menhaden_Catch_StationarySystem.csv",
"F.menhaden_B_StationarySystem.csv"
),
reference_points_scenario = "compensation",
functional_groups = functional_groups,
key_functional_group = "AtlanticMenhaden",
ages = ages,
maturity_vec = maturity_matrix[1,]
)
stationary_fmsy_s1_CatchVec <- stationary_msy_s1_MaxCatchVec$FMSY * stationary_fmsy_s1[[1]]$Fbase
msy_file_path <- here::here("data", "ewe", "7ages_ecosim_om_msy", "msy_forcing_pdsi_egg_amo1_fleet_details")
compensation_msy_s2_MaxSumCatch <- 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,
maturity_vec = maturity_matrix[1,]
)
compensation_fmsy_s2_MaxSumCatch <- compensation_msy_s2_MaxSumCatch$FMSY * compensation_fmsy_s2[[1]]$Fbase
compensation_msy_s2_MaxCatchVec <- read_ewe_reference_points_vec(
file_path = msy_file_path,
file_names = c(
"F.menhaden_Catch_FullCompensation.csv",
"F.menhaden_B_FullCompensation.csv"
),
reference_points_scenario = "compensation",
functional_groups = functional_groups,
key_functional_group = "AtlanticMenhaden",
ages = ages,
maturity_vec = maturity_matrix[1,]
)
compensation_fmsy_s2_CatchVec <- compensation_msy_s2_MaxCatchVec$FMSY * compensation_fmsy_s2[[1]]$Fbase
stationary_msy_s2_MaxSumCatch <- 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,
maturity_vec = maturity_matrix[1,]
)
stationary_fmsy_s2_MaxSumCatch <- stationary_msy_s2_MaxSumCatch$FMSY * stationary_fmsy_s2[[1]]$Fbase
stationary_msy_s2_MaxCatchVec <- read_ewe_reference_points_vec(
file_path = msy_file_path,
file_names = c(
"F.menhaden_Catch_StationarySystem.csv",
"F.menhaden_B_StationarySystem.csv"
),
reference_points_scenario = "compensation",
functional_groups = functional_groups,
key_functional_group = "AtlanticMenhaden",
ages = ages,
maturity_vec = maturity_matrix[1,]
)
stationary_fmsy_s2_CatchVec <- stationary_msy_s2_MaxCatchVec$FMSY * stationary_fmsy_s2[[1]]$Fbase
fmsy_table <- cbind( # compensation_fmsy_s1_MaxSumCatch,
compensation_fmsy_s1_CatchVec,
# stationary_fmsy_s1_MaxSumCatch,
stationary_fmsy_s1_CatchVec,
# compensation_fmsy_s2_MaxSumCatch,
compensation_fmsy_s2_CatchVec,
# stationary_fmsy_s2_MaxSumCatch,
stationary_fmsy_s2_CatchVec
)
fmsy_table <- rbind(fmsy_table, apply(fmsy_table, 2, max))
row.names(fmsy_table) <- c(paste("Case stock", c(0:5, "6+")), "Apical FMSY")
knitr::kable(fmsy_table,
caption = "Table 5. FMSY output from catch over F output"
)
msy_table <- c(
# compensation_msy_s1_MaxSumCatch$MSY *1000000,
sum(compensation_msy_s1_MaxCatchVec$MSY) * 1000000,
# stationary_msy_s1_MaxSumCatch$MSY *1000000,
sum(stationary_msy_s1_MaxCatchVec$MSY) * 1000000,
# compensation_msy_s2_MaxSumCatch$MSY *1000000,
sum(compensation_msy_s2_MaxCatchVec$MSY) * 1000000,
# stationary_msy_s2_MaxSumCatch$MSY*1000000,
sum(stationary_msy_s2_MaxCatchVec$MSY) * 1000000
)
names(msy_table) <- c(
# "compensation_msy_s1_MaxSumCatch",
"Compensation MSY s1",
# "stationary_msy_s1_MaxSumCatch",
"Stationary MSY s1",
# "compensation_msy_s2_MaxSumCatch",
"Compensation msy s2",
# "stationary_msy_s2_MaxSumCatch",
"Stationary MSY s2"
)
knitr::kable(msy_table,
caption = "Table 5. MSY output from catch over F output"
)
bmsy_table <- c(
# compensation_msy_s1_MaxSumCatch$BMSY *1000000,
sum(compensation_msy_s1_MaxCatchVec$BMSY) * 1000000,
# stationary_msy_s1_MaxSumCatch$BMSY *1000000,
sum(stationary_msy_s1_MaxCatchVec$BMSY) * 1000000,
# compensation_msy_s2_MaxSumCatch$BMSY *1000000,
sum(compensation_msy_s2_MaxCatchVec$BMSY) * 1000000,
# stationary_msy_s2_MaxSumCatch$BMSY*1000000,
sum(stationary_msy_s2_MaxCatchVec$BMSY) * 1000000
)
names(bmsy_table) <- c(
# "compensation_bmsy_s1_MaxSumCatch",
"Compensation BMSY s1",
# "stationary_bmsy_s1_MaxSumCatch",
"Stationary BMSY s1",
# "compensation_bmsy_s2_MaxSumCatch",
"Compensation BMSY s2",
# "stationary_bmsy_s2_MaxSumCatch",
"Stationary BMSY s2"
)
sbmsy_table <- c(
sum(compensation_msy_s1_MaxCatchVec$SBMSY) * 1000000,
sum(stationary_msy_s1_MaxCatchVec$SBMSY) * 1000000,
sum(compensation_msy_s2_MaxCatchVec$SBMSY) * 1000000,
sum(stationary_msy_s2_MaxCatchVec$SBMSY) * 1000000
)
names(sbmsy_table) <- c(
"Compensation SBMSY s1",
"Stationary SBMSY s1",
"Compensation SBMSY s2",
"Stationary SBMSY s2"
)
knitr::kable(sbmsy_table,
caption = "Table 6. SBMSY output from catch over F output"
)
# Bratio and SBratio
biomass_s1 <- read_ewe_output(
file_path = here::here(
"data", "ewe", "7ages_ecosim_om",
"ecosim_base_run"
),
file_names = "biomass_monthly.csv",
skip_nrows = 8,
plot = FALSE,
figure_titles = NULL,
functional_groups = functional_groups,
figure_colors = NULL
)
model_year <- 1985:2012
time_id <- seq(1, nrow(biomass_s1[[1]]), by = 12)[1:(length(model_year))]
biomass_ewe_s1 <- apply(biomass_s1[[1]][, age_name], 1, sum) * 1000000
sb_ewe_s1 <- apply(as.matrix(biomass_s1[[1]][, age_name]) %*% diag(maturity_matrix[1, ]), 1, sum) * 1000000
biomass_s2 <- read_ewe_output(
file_path = here::here(
"data", "ewe", "7ages_ecosim_om",
"ecosim_forcing_pdsi_egg_amo1"
),
file_names = "biomass_monthly.csv",
skip_nrows = 8,
plot = FALSE,
figure_titles = NULL,
functional_groups = functional_groups,
figure_colors = NULL
)
biomass_ewe_s2 <- apply(biomass_s2[[1]][, age_name], 1, sum) * 1000000
sb_ewe_s2 <- apply(as.matrix(biomass_s2[[1]][, age_name]) %*% diag(maturity_matrix[1, ]), 1, sum) * 1000000
biomass_s3 <- read_ewe_output(
file_path = here::here(
"data", "ewe", "7ages_ecosim_om",
"ecosim_fleet_dynamics"
),
file_names = "biomass_monthly.csv",
skip_nrows = 8,
plot = FALSE,
figure_titles = NULL,
functional_groups = functional_groups,
figure_colors = NULL
)
biomass_ewe_s3 <- apply(biomass_s3[[1]][, age_name], 1, sum) * 1000000
sb_ewe_s3 <- apply(as.matrix(biomass_s3[[1]][, age_name]) %*% diag(maturity_matrix[1, ]), 1, sum) * 1000000
bratio_table <- c(
# biomass_ewe_s1[tail(time_id, n=1)] / (compensation_msy_s1_MaxSumCatch$BMSY *1000000),
biomass_ewe_s1[tail(time_id, n = 1)] / (sum(compensation_msy_s1_MaxCatchVec$BMSY) * 1000000),
# biomass_ewe_s1[tail(time_id, n=1)] / (stationary_msy_s1_MaxSumCatch$BMSY *1000000),
biomass_ewe_s1[tail(time_id, n = 1)] / (sum(stationary_msy_s1_MaxCatchVec$BMSY) * 1000000),
# biomass_ewe_s2[tail(time_id, n=1)] / (compensation_msy_s2_MaxSumCatch$BMSY *1000000),
biomass_ewe_s2[tail(time_id, n = 1)] / (sum(compensation_msy_s2_MaxCatchVec$BMSY) * 1000000),
# biomass_ewe_s2[tail(time_id, n=1)] / (stationary_msy_s2_MaxSumCatch$BMSY*1000000),
biomass_ewe_s2[tail(time_id, n = 1)] / (sum(stationary_msy_s2_MaxCatchVec$BMSY) * 1000000),
# biomass_ewe_s3[tail(time_id, n=1)] / (compensation_msy_s2_MaxSumCatch$BMSY *1000000),
biomass_ewe_s3[tail(time_id, n = 1)] / (sum(compensation_msy_s2_MaxCatchVec$BMSY) * 1000000),
# biomass_ewe_s3[tail(time_id, n=1)] / (stationary_msy_s2_MaxSumCatch$BMSY*1000000),
biomass_ewe_s3[tail(time_id, n = 1)] / (sum(stationary_msy_s2_MaxCatchVec$BMSY) * 1000000)
)
names(bratio_table) <- c(
# "compensation_bratio_s1_MaxSumCatch",
"compensation_bratio_s1",
# "stationary_bratio_s1_MaxSumCatch",
"stationary_bratio_s1",
# "compensation_bratio_s2_MaxSumCatch",
"compensation_bratio_s2",
# "stationary_bratio_s2_MaxSumCatch",
"stationary_bratio_s2",
# "compensation_bratio_s3_MaxSumCatch",
"compensation_bratio_s3",
# "stationary_bratio_s3_MaxSumCatch",
"stationary_bratio_s3"
)
knitr::kable(bratio_table,
caption = "Table 7. Bratio output from catch over F output"
)
sbratio_table <- c(
sb_ewe_s1[tail(time_id, n = 1)] / (sum(compensation_msy_s1_MaxCatchVec$SBMSY) * 1000000),
sb_ewe_s1[tail(time_id, n = 1)] / (sum(stationary_msy_s1_MaxCatchVec$SBMSY) * 1000000),
sb_ewe_s2[tail(time_id, n = 1)] / (sum(compensation_msy_s2_MaxCatchVec$SBMSY) * 1000000),
sb_ewe_s2[tail(time_id, n = 1)] / (sum(stationary_msy_s2_MaxCatchVec$SBMSY) * 1000000),
sb_ewe_s3[tail(time_id, n = 1)] / (sum(compensation_msy_s2_MaxCatchVec$SBMSY) * 1000000),
sb_ewe_s3[tail(time_id, n = 1)] / (sum(stationary_msy_s2_MaxCatchVec$SBMSY) * 1000000)
)
names(sbratio_table) <- c(
"compensation_sbratio_s1",
"stationary_sbratio_s1",
"compensation_sbratio_s2",
"stationary_sbratio_s2",
"compensation_sbratio_s3",
"stationary_sbratio_s3"
)
knitr::kable(sbratio_table,
caption = "Table 8. SBratio output from catch over F output"
)
reference_points <- list(
fmsy=fmsy_table,
bmsy=bmsy_table,
msy=msy_table,
bratio=bratio_table,
sbratio=sbratio_table
)
save(reference_points, file = here::here("data", "ewe", "ewe_reference_points.RData"))
Bai-Li-NOAA/IFA4EBFM documentation built on July 1, 2024, 4:53 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( echo = FALSE, message = FALSE, warning = FALSE, fig.width = 12, fig.height = 10, results = "asis" ) options(scipen = 1, digits = 4)
devtools::load_all() library(tidyr)
Evaluation of equilibrium MSY reference points from EwE user manual
- Single-species FMSY (stationary FMSY)
- Run the Ecosim model to equilibrium for a range of F values while holding biomasses of all other groups constant (except other life history stanzas for the same species in split and multi-stanza cases)
- Examine predicted compensatory responses by the group caused by the foraging arena functional response and related foraging time adjustment parameters
- Ecological FMSY (compensatory FMSY)
- Allows all groups' biomass to change in response to the biomass changes of the focal FMSY groups
- The assessment usually produces a more optimistic estimate of MSY for each group
- Increases in food species abundances and decreases in predator abundances when any one group is harvested harder
Output comparison
Run model to get compensation FMSY vs stationary FMSY
# load sa_data to get maturity data scenario_filename <- "_lowCV" terminal_year <- 2012 fishery_CV_input <- 0.05 survey_CV_input <- 0.1 add_environmental_effects <- FALSE add_fleet_dynamics <- FALSE source(here::here("Rscript", "simulationOM3.R")) maturity_matrix <- sa_data$biodata$maturity_matrix # load 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) # base scenario compensation_fmsy_s1 <- read_fmsy_table( file_names = here::here( "data", "ewe", "7ages_ecosim_om_msy", paste0("msy_base_run_age", ages), "FMSY_fullCompensation.csv" ), colnames = "Fmsy", functional_groups = functional_groups ) compensation_fmsy_s1 <- c( compensation_fmsy_s1, read_fmsy_table( file_names = here::here( "data", "ewe", "7ages_ecosim_om_msy", "msy_base_run_fleet", "FMSY_fullCompensation.csv" ), colnames = "Fmsy", functional_groups = functional_groups ) ) stationary_fmsy_s1 <- read_fmsy_table( file_names = here::here( "data", "ewe", "7ages_ecosim_om_msy", paste0("msy_base_run_age", ages), "FMSY_StationarySystem.csv" ), colnames = "Fmsy", functional_groups = functional_groups ) stationary_fmsy_s1 <- c( stationary_fmsy_s1, read_fmsy_table( file_names = here::here( "data", "ewe", "7ages_ecosim_om_msy", "msy_base_run_fleet", "FMSY_StationarySystem.csv" ), colnames = "Fmsy", functional_groups = functional_groups ) ) mean_difference <- all.equal(current = stationary_fmsy_s1, target = compensation_fmsy_s1) mean_difference_value <- as.numeric(sapply(strsplit(mean_difference, " "), "[[", 6)) var_s1 <- var(mean_difference_value) # Environmental forcing scenario compensation_fmsy_s2 <- read_fmsy_table( file_names = here::here( "data", "ewe", "7ages_ecosim_om_msy", paste0("msy_forcing_pdsi_egg_amo1_age", ages), "FMSY_fullCompensation.csv" ), colnames = "Fmsy", functional_groups = functional_groups ) compensation_fmsy_s2 <- c( compensation_fmsy_s2, read_fmsy_table( file_names = here::here( "data", "ewe", "7ages_ecosim_om_msy", "msy_forcing_pdsi_egg_amo1_fleet", "FMSY_fullCompensation.csv" ), colnames = "Fmsy", functional_groups = functional_groups ) ) stationary_fmsy_s2 <- read_fmsy_table( file_names = here::here( "data", "ewe", "7ages_ecosim_om_msy", paste0("msy_forcing_pdsi_egg_amo1_age", ages), "FMSY_StationarySystem.csv" ), colnames = "Fmsy", functional_groups = functional_groups ) stationary_fmsy_s2 <- c( stationary_fmsy_s2, read_fmsy_table( file_names = here::here( "data", "ewe", "7ages_ecosim_om_msy", "msy_forcing_pdsi_egg_amo1_fleet", "FMSY_StationarySystem.csv" ), colnames = "Fmsy", functional_groups = functional_groups ) ) mean_difference <- all.equal(current = stationary_fmsy_s2, target = compensation_fmsy_s2) mean_difference_value <- as.numeric(sapply(strsplit(mean_difference, " "), "[[", 6)) var_s2 <- var(mean_difference_value) if (var_s1 == 0 & var_s2 == 0) { cat("- Mean relative differences between compensation and stationary FMSYs across 7 age classes from scenarios S1 and S2 = ", unique(mean_difference_value), "\n") } else { cat("- Mean relative differences between compensation and stationary FMSYs vary by age groups\n") }
Run model by functional group vs by fleet
mean_difference <- c() list_length <- length(compensation_fmsy_s1) # base scenario for (i in seq_along(compensation_fmsy_s1)) { mean_difference[i] <- all.equal(current = compensation_fmsy_s1[[list_length]], target = compensation_fmsy_s1[[i]]) } compensation_mean_difference_s1 <- all(mean_difference) for (i in seq_along(stationary_fmsy_s1)) { mean_difference[i] <- all.equal(current = stationary_fmsy_s1[[list_length]], target = stationary_fmsy_s1[[i]]) } stationary_mean_difference_s1 <- all(mean_difference) # environmental forcing scenario for (i in seq_along(compensation_fmsy_s2)) { mean_difference[i] <- all.equal(current = compensation_fmsy_s2[[list_length]], target = compensation_fmsy_s2[[i]]) } compensation_mean_difference_s2 <- all(mean_difference) for (i in seq_along(stationary_fmsy_s2)) { mean_difference[i] <- all.equal(current = stationary_fmsy_s2[[list_length]], target = stationary_fmsy_s2[[i]]) } stationary_mean_difference_s2 <- all(mean_difference) if (all(c(compensation_mean_difference_s1, stationary_mean_difference_s1, compensation_mean_difference_s2, stationary_mean_difference_s2))) { cat("- FMSYs by functional group are identical to FMSYs by fleet\n") } else { cat("- FMSYs by functional group are different to FMSYs by fleet\n") }
Run model to equilibrium for a range of F values vs run to depletion
msy_ages_label <- c(0:5, "6+") maxF_file_path <- here::here( "data", "ewe", "7ages_ecosim_om_msy", paste0("msy_forcing_pdsi_egg_amo1_age", ages, "_details") ) maxF_file_path <- c( maxF_file_path, here::here("data", "ewe", "7ages_ecosim_om_msy", "msy_forcing_pdsi_egg_amo1_fleet_details") ) catch_file_names <- c(paste0("menhaden ", msy_ages_label, "_Catch_FullCompensation.csv"), "F.menhaden_Catch_FullCompensation.csv") biomass_file_names <- c(paste0("menhaden ", msy_ages_label, "_B_FullCompensation.csv"), "F.menhaden_B_FullCompensation.csv") compensation_msy_maxF <- list() for (i in seq_along(maxF_file_path)) { compensation_msy_maxF[[i]] <- read_ewe_reference_points( file_path = maxF_file_path[i], file_names = c( catch_file_names[i], biomass_file_names[i] ), reference_points_scenario = "compensation", functional_groups = functional_groups, key_functional_group = "AtlanticMenhaden", ages = ages, plot = FALSE, maturity_vec = maturity_matrix[1,] ) } depletion_file_path <- here::here( "data", "ewe", "7ages_ecosim_om_msy", paste0("msy_forcing_pdsi_egg_amo1_age", ages, "_details_depletion") ) depletion_file_path <- c( depletion_file_path, here::here("data", "ewe", "7ages_ecosim_om_msy", "msy_forcing_pdsi_egg_amo1_fleet_details_depletion") ) compensation_msy_depletion <- list() for (i in seq_along(maxF_file_path)) { compensation_msy_depletion[[i]] <- read_ewe_reference_points( file_path = depletion_file_path[i], file_names = c( catch_file_names[i], biomass_file_names[i] ), reference_points_scenario = "compensation", functional_groups = functional_groups, key_functional_group = "AtlanticMenhaden", ages = ages, plot = FALSE, maturity_vec = maturity_matrix[1,] ) } result_compensation <- all.equal(current = compensation_msy_maxF, target = compensation_msy_depletion) catch_file_names <- c(paste0("menhaden ", msy_ages_label, "_Catch_StationarySystem.csv"), "F.menhaden_Catch_StationarySystem.csv") biomass_file_names <- c(paste0("menhaden ", msy_ages_label, "_B_StationarySystem.csv"), "F.menhaden_B_StationarySystem.csv") stationary_msy_maxF <- list() for (i in seq_along(maxF_file_path)) { stationary_msy_maxF[[i]] <- read_ewe_reference_points( file_path = maxF_file_path[i], file_names = c( catch_file_names[i], biomass_file_names[i] ), reference_points_scenario = "stationary", functional_groups = functional_groups, key_functional_group = "AtlanticMenhaden", ages = ages, plot = FALSE, maturity_vec = maturity_matrix[1,] ) } stationary_msy_depletion <- list() for (i in seq_along(maxF_file_path)) { stationary_msy_depletion[[i]] <- read_ewe_reference_points( file_path = depletion_file_path[i], file_names = c( catch_file_names[i], biomass_file_names[i] ), reference_points_scenario = "stationary", functional_groups = functional_groups, key_functional_group = "AtlanticMenhaden", ages = ages, plot = FALSE, maturity_vec = maturity_matrix[1,] ) } result_stationary <- all.equal(current = stationary_msy_maxF, target = stationary_msy_depletion) if (all(c(result_compensation, result_stationary))) { cat("- FMSYs from the two searching methods are identical\n") } else { cat("- FMSYs from the two searching methods are different\n") }
Single FMSY vs multiple FMSYs
- EwE default summary output files (e.g., FMSY_FullCompensation.csv or FMSY_StationarySystem.csv) show the estimate of FMSY for each functional group (Table 1-4).
- EwE also produces the catch of each functional group over a range of F values in F.menhaden_Catch_FullCompensation.csv or F.menhaden_Catch_StationarySystem.csv. To extract reference points, we can get a vector of fishing effort multipliers that correspond to the maximum catch of the case stock, and them multiply the fishing effort multipliers with seven base fishing mortality rates from 1985 to obtain seven FMSY values by age.
compensation_fmsy_s1 <- read_fmsy_table( file_names = here::here( "data", "ewe", "7ages_ecosim_om_msy", "msy_base_run_fleet", "FMSY_fullCompensation.csv" ), colnames = c( "Group", "TL", "Fbase", "Cbase", "Vbase", "FmsyFound", "Fmsy", "Cmsy", "Vmsy" ), functional_groups = functional_groups ) stationary_fmsy_s1 <- read_fmsy_table( file_names = here::here( "data", "ewe", "7ages_ecosim_om_msy", "msy_base_run_fleet", "FMSY_StationarySystem.csv" ), colnames = c( "Group", "TL", "Fbase", "Cbase", "Vbase", "FmsyFound", "Fmsy", "Cmsy", "Vmsy" ), functional_groups = functional_groups ) knitr::kable(compensation_fmsy_s1, caption = "Table 1. Base scenario S1 compensation FMSY output from EwE" ) knitr::kable(stationary_fmsy_s1, caption = "Table 2. Base scenario S1 stationary FMSY output from EwE" ) compensation_fmsy_s2 <- read_fmsy_table( file_names = here::here( "data", "ewe", "7ages_ecosim_om_msy", "msy_forcing_pdsi_egg_amo1_fleet", "FMSY_fullCompensation.csv" ), colnames = c( "Group", "TL", "Fbase", "Cbase", "Vbase", "FmsyFound", "Fmsy", "Cmsy", "Vmsy" ), functional_groups = functional_groups ) stationary_fmsy_s2 <- read_fmsy_table( file_names = here::here( "data", "ewe", "7ages_ecosim_om_msy", "msy_forcing_pdsi_egg_amo1_fleet", "FMSY_StationarySystem.csv" ), colnames = c( "Group", "TL", "Fbase", "Cbase", "Vbase", "FmsyFound", "Fmsy", "Cmsy", "Vmsy" ), functional_groups = functional_groups ) knitr::kable(compensation_fmsy_s2, caption = "Table 3. Environmental forcing scenario S2 compensation FMSY output from EwE" ) knitr::kable(stationary_fmsy_s2, caption = "Table 4. Environmental forcing scenario S2 stationary FMSY output from EwE" ) msy_file_path <- here::here("data", "ewe", "7ages_ecosim_om_msy", "msy_base_run_fleet_details") compensation_msy_s1_MaxSumCatch <- 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, maturity_vec = maturity_matrix[1,] ) compensation_fmsy_s1_MaxSumCatch <- compensation_msy_s1_MaxSumCatch$FMSY * compensation_fmsy_s1[[1]]$Fbase compensation_msy_s1_MaxCatchVec <- read_ewe_reference_points_vec( file_path = msy_file_path, file_names = c( "F.menhaden_Catch_FullCompensation.csv", "F.menhaden_B_FullCompensation.csv" ), reference_points_scenario = "compensation", functional_groups = functional_groups, key_functional_group = "AtlanticMenhaden", ages = ages, maturity_vec = maturity_matrix[1,] ) compensation_fmsy_s1_CatchVec <- compensation_msy_s1_MaxCatchVec$FMSY * compensation_fmsy_s1[[1]]$Fbase stationary_msy_s1_MaxSumCatch <- 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, maturity_vec = maturity_matrix[1,] ) stationary_fmsy_s1_MaxSumCatch <- stationary_msy_s1_MaxSumCatch$FMSY * stationary_fmsy_s1[[1]]$Fbase stationary_msy_s1_MaxCatchVec <- read_ewe_reference_points_vec( file_path = msy_file_path, file_names = c( "F.menhaden_Catch_StationarySystem.csv", "F.menhaden_B_StationarySystem.csv" ), reference_points_scenario = "compensation", functional_groups = functional_groups, key_functional_group = "AtlanticMenhaden", ages = ages, maturity_vec = maturity_matrix[1,] ) stationary_fmsy_s1_CatchVec <- stationary_msy_s1_MaxCatchVec$FMSY * stationary_fmsy_s1[[1]]$Fbase msy_file_path <- here::here("data", "ewe", "7ages_ecosim_om_msy", "msy_forcing_pdsi_egg_amo1_fleet_details") compensation_msy_s2_MaxSumCatch <- 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, maturity_vec = maturity_matrix[1,] ) compensation_fmsy_s2_MaxSumCatch <- compensation_msy_s2_MaxSumCatch$FMSY * compensation_fmsy_s2[[1]]$Fbase compensation_msy_s2_MaxCatchVec <- read_ewe_reference_points_vec( file_path = msy_file_path, file_names = c( "F.menhaden_Catch_FullCompensation.csv", "F.menhaden_B_FullCompensation.csv" ), reference_points_scenario = "compensation", functional_groups = functional_groups, key_functional_group = "AtlanticMenhaden", ages = ages, maturity_vec = maturity_matrix[1,] ) compensation_fmsy_s2_CatchVec <- compensation_msy_s2_MaxCatchVec$FMSY * compensation_fmsy_s2[[1]]$Fbase stationary_msy_s2_MaxSumCatch <- 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, maturity_vec = maturity_matrix[1,] ) stationary_fmsy_s2_MaxSumCatch <- stationary_msy_s2_MaxSumCatch$FMSY * stationary_fmsy_s2[[1]]$Fbase stationary_msy_s2_MaxCatchVec <- read_ewe_reference_points_vec( file_path = msy_file_path, file_names = c( "F.menhaden_Catch_StationarySystem.csv", "F.menhaden_B_StationarySystem.csv" ), reference_points_scenario = "compensation", functional_groups = functional_groups, key_functional_group = "AtlanticMenhaden", ages = ages, maturity_vec = maturity_matrix[1,] ) stationary_fmsy_s2_CatchVec <- stationary_msy_s2_MaxCatchVec$FMSY * stationary_fmsy_s2[[1]]$Fbase fmsy_table <- cbind( # compensation_fmsy_s1_MaxSumCatch, compensation_fmsy_s1_CatchVec, # stationary_fmsy_s1_MaxSumCatch, stationary_fmsy_s1_CatchVec, # compensation_fmsy_s2_MaxSumCatch, compensation_fmsy_s2_CatchVec, # stationary_fmsy_s2_MaxSumCatch, stationary_fmsy_s2_CatchVec ) fmsy_table <- rbind(fmsy_table, apply(fmsy_table, 2, max)) row.names(fmsy_table) <- c(paste("Case stock", c(0:5, "6+")), "Apical FMSY") knitr::kable(fmsy_table, caption = "Table 5. FMSY output from catch over F output" ) msy_table <- c( # compensation_msy_s1_MaxSumCatch$MSY *1000000, sum(compensation_msy_s1_MaxCatchVec$MSY) * 1000000, # stationary_msy_s1_MaxSumCatch$MSY *1000000, sum(stationary_msy_s1_MaxCatchVec$MSY) * 1000000, # compensation_msy_s2_MaxSumCatch$MSY *1000000, sum(compensation_msy_s2_MaxCatchVec$MSY) * 1000000, # stationary_msy_s2_MaxSumCatch$MSY*1000000, sum(stationary_msy_s2_MaxCatchVec$MSY) * 1000000 ) names(msy_table) <- c( # "compensation_msy_s1_MaxSumCatch", "Compensation MSY s1", # "stationary_msy_s1_MaxSumCatch", "Stationary MSY s1", # "compensation_msy_s2_MaxSumCatch", "Compensation msy s2", # "stationary_msy_s2_MaxSumCatch", "Stationary MSY s2" ) knitr::kable(msy_table, caption = "Table 5. MSY output from catch over F output" ) bmsy_table <- c( # compensation_msy_s1_MaxSumCatch$BMSY *1000000, sum(compensation_msy_s1_MaxCatchVec$BMSY) * 1000000, # stationary_msy_s1_MaxSumCatch$BMSY *1000000, sum(stationary_msy_s1_MaxCatchVec$BMSY) * 1000000, # compensation_msy_s2_MaxSumCatch$BMSY *1000000, sum(compensation_msy_s2_MaxCatchVec$BMSY) * 1000000, # stationary_msy_s2_MaxSumCatch$BMSY*1000000, sum(stationary_msy_s2_MaxCatchVec$BMSY) * 1000000 ) names(bmsy_table) <- c( # "compensation_bmsy_s1_MaxSumCatch", "Compensation BMSY s1", # "stationary_bmsy_s1_MaxSumCatch", "Stationary BMSY s1", # "compensation_bmsy_s2_MaxSumCatch", "Compensation BMSY s2", # "stationary_bmsy_s2_MaxSumCatch", "Stationary BMSY s2" ) sbmsy_table <- c( sum(compensation_msy_s1_MaxCatchVec$SBMSY) * 1000000, sum(stationary_msy_s1_MaxCatchVec$SBMSY) * 1000000, sum(compensation_msy_s2_MaxCatchVec$SBMSY) * 1000000, sum(stationary_msy_s2_MaxCatchVec$SBMSY) * 1000000 ) names(sbmsy_table) <- c( "Compensation SBMSY s1", "Stationary SBMSY s1", "Compensation SBMSY s2", "Stationary SBMSY s2" ) knitr::kable(sbmsy_table, caption = "Table 6. SBMSY output from catch over F output" ) # Bratio and SBratio biomass_s1 <- read_ewe_output( file_path = here::here( "data", "ewe", "7ages_ecosim_om", "ecosim_base_run" ), file_names = "biomass_monthly.csv", skip_nrows = 8, plot = FALSE, figure_titles = NULL, functional_groups = functional_groups, figure_colors = NULL ) model_year <- 1985:2012 time_id <- seq(1, nrow(biomass_s1[[1]]), by = 12)[1:(length(model_year))] biomass_ewe_s1 <- apply(biomass_s1[[1]][, age_name], 1, sum) * 1000000 sb_ewe_s1 <- apply(as.matrix(biomass_s1[[1]][, age_name]) %*% diag(maturity_matrix[1, ]), 1, sum) * 1000000 biomass_s2 <- read_ewe_output( file_path = here::here( "data", "ewe", "7ages_ecosim_om", "ecosim_forcing_pdsi_egg_amo1" ), file_names = "biomass_monthly.csv", skip_nrows = 8, plot = FALSE, figure_titles = NULL, functional_groups = functional_groups, figure_colors = NULL ) biomass_ewe_s2 <- apply(biomass_s2[[1]][, age_name], 1, sum) * 1000000 sb_ewe_s2 <- apply(as.matrix(biomass_s2[[1]][, age_name]) %*% diag(maturity_matrix[1, ]), 1, sum) * 1000000 biomass_s3 <- read_ewe_output( file_path = here::here( "data", "ewe", "7ages_ecosim_om", "ecosim_fleet_dynamics" ), file_names = "biomass_monthly.csv", skip_nrows = 8, plot = FALSE, figure_titles = NULL, functional_groups = functional_groups, figure_colors = NULL ) biomass_ewe_s3 <- apply(biomass_s3[[1]][, age_name], 1, sum) * 1000000 sb_ewe_s3 <- apply(as.matrix(biomass_s3[[1]][, age_name]) %*% diag(maturity_matrix[1, ]), 1, sum) * 1000000 bratio_table <- c( # biomass_ewe_s1[tail(time_id, n=1)] / (compensation_msy_s1_MaxSumCatch$BMSY *1000000), biomass_ewe_s1[tail(time_id, n = 1)] / (sum(compensation_msy_s1_MaxCatchVec$BMSY) * 1000000), # biomass_ewe_s1[tail(time_id, n=1)] / (stationary_msy_s1_MaxSumCatch$BMSY *1000000), biomass_ewe_s1[tail(time_id, n = 1)] / (sum(stationary_msy_s1_MaxCatchVec$BMSY) * 1000000), # biomass_ewe_s2[tail(time_id, n=1)] / (compensation_msy_s2_MaxSumCatch$BMSY *1000000), biomass_ewe_s2[tail(time_id, n = 1)] / (sum(compensation_msy_s2_MaxCatchVec$BMSY) * 1000000), # biomass_ewe_s2[tail(time_id, n=1)] / (stationary_msy_s2_MaxSumCatch$BMSY*1000000), biomass_ewe_s2[tail(time_id, n = 1)] / (sum(stationary_msy_s2_MaxCatchVec$BMSY) * 1000000), # biomass_ewe_s3[tail(time_id, n=1)] / (compensation_msy_s2_MaxSumCatch$BMSY *1000000), biomass_ewe_s3[tail(time_id, n = 1)] / (sum(compensation_msy_s2_MaxCatchVec$BMSY) * 1000000), # biomass_ewe_s3[tail(time_id, n=1)] / (stationary_msy_s2_MaxSumCatch$BMSY*1000000), biomass_ewe_s3[tail(time_id, n = 1)] / (sum(stationary_msy_s2_MaxCatchVec$BMSY) * 1000000) ) names(bratio_table) <- c( # "compensation_bratio_s1_MaxSumCatch", "compensation_bratio_s1", # "stationary_bratio_s1_MaxSumCatch", "stationary_bratio_s1", # "compensation_bratio_s2_MaxSumCatch", "compensation_bratio_s2", # "stationary_bratio_s2_MaxSumCatch", "stationary_bratio_s2", # "compensation_bratio_s3_MaxSumCatch", "compensation_bratio_s3", # "stationary_bratio_s3_MaxSumCatch", "stationary_bratio_s3" ) knitr::kable(bratio_table, caption = "Table 7. Bratio output from catch over F output" ) sbratio_table <- c( sb_ewe_s1[tail(time_id, n = 1)] / (sum(compensation_msy_s1_MaxCatchVec$SBMSY) * 1000000), sb_ewe_s1[tail(time_id, n = 1)] / (sum(stationary_msy_s1_MaxCatchVec$SBMSY) * 1000000), sb_ewe_s2[tail(time_id, n = 1)] / (sum(compensation_msy_s2_MaxCatchVec$SBMSY) * 1000000), sb_ewe_s2[tail(time_id, n = 1)] / (sum(stationary_msy_s2_MaxCatchVec$SBMSY) * 1000000), sb_ewe_s3[tail(time_id, n = 1)] / (sum(compensation_msy_s2_MaxCatchVec$SBMSY) * 1000000), sb_ewe_s3[tail(time_id, n = 1)] / (sum(stationary_msy_s2_MaxCatchVec$SBMSY) * 1000000) ) names(sbratio_table) <- c( "compensation_sbratio_s1", "stationary_sbratio_s1", "compensation_sbratio_s2", "stationary_sbratio_s2", "compensation_sbratio_s3", "stationary_sbratio_s3" ) knitr::kable(sbratio_table, caption = "Table 8. SBratio output from catch over F output" ) reference_points <- list( fmsy=fmsy_table, bmsy=bmsy_table, msy=msy_table, bratio=bratio_table, sbratio=sbratio_table ) save(reference_points, file = here::here("data", "ewe", "ewe_reference_points.RData"))
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.