In Bai-Li-NOAA/IFA4EBFM: Interdisciplinary Projections for Ecosystem-Based Fisheries Management
knitr::opts_chunk$set(echo = TRUE)
The Ecopath with Ecosim (EwE) Model
The goal of the this work is to restructure an EwE model that include a menhaden-like species with 7 age classes (i.e. 0 - 6+). An existing EwE model of the Northwest Atlantic Continental Shelf model of intermediate complexity for ecosystems (NWACS-MICE) was previously developed to link the dynamics of menhaden with key managed predators (Chagaris et al., 2020). There were two age classes of menhaden in the NWACS-MICE model. To restructure the MWACS-MICE model, we made following changes:
-
Updated Ecopath multi-stanza groups (menhaden-like species)
-
Updated age 6+ biomass in 1985 using Beaufort Assessment Model (BAM) outputs (SEDAR 2020).
-
Updated total mortality-at-age in 1985 inputs using $C_{a}/B_{a}+M_{a}$ with BAM outputs. Here $C$, $B$, and $M$ represent catch, biomass, and natural mortality respectively. $a$ denotes age.
-
Updated relative biomass accumulation rate (BA/B). Three scenarios were explored: 1) 0.114 from Chagaris et al. (2020), 2) EwE default value of 0, 3) -0.081 based on $(B_{1986} - B_{1985})/B_{1986}$ (Buchheister et al., 2017). Both BA/B = 0 and BA/B = -0.081 help improve agreement between Ecosim predicted catch with observed catch and BA/B = 0 showed lower sum of squares (SS) value. SS values from scenarios 1 - 3 were 3059, 2339, and 3136 respectively. So we used BA/B = 0 for base run. Key figures from the three scenarios could be found here.
-
Updated Ecopath diet composition using inputs from NWACS-MICE model (Chagaris et al., 2020), full NWACS model (Buchhester et al., 2017), and BAM M-at-age data.
-
Updated Ecopath C-at-age using BAM outputs.
-
Updated Ecosim vulnerabilities and ages 1-6+ share the same vulnerability values.
-
Updated time series matrix using I-at-age, C-at-age, and F-at-age from BAM outputs.
Here $I$ and $F$ represent abundance indices and fishing mortality respectively.
-
A few forcing scenarios were tested to explore various sources of environmental variability. For example, primary production splines was applied to phytoplankton, recruitment deviations from Chagaris et al. (2020) and BAM outputs were linked with menhaden age 0 fish. However, incorporating primary production splines did help soak up some of the process variability in the food web dynamics and did not help improve model fits. Incorporation of recruitment deviations did not help better represent interannual variability in the system. Question: Maybe not apply forcing functions to the EwE model, but exploring relationships between environmental drivers and menhaden biomass when adjusting $F_{ECO}$?.
-
Density-dependent growth: According to EwE 6.0 guide, users could set a non-zero feeding time adjustment for the juvenile group and this represents density-dependent changes in juvenile mortality rate associated with changes in feeding time and predation risk. The feeding time adjust rate for age 0 of menhaden-like species was set to 0.5.
Comparison between EwE model with 2 age classes of menhaden (Chagaris et al., 2020) and EwE model with 7 age classes of menhaden-like species.
devtools::load_all()
# Load output data from BAM and EwE ---------------------------------------
project_path <- here::here()
ewe2ages_path <- file.path(project_path, "data", "ewe", "2ages")
ewe7ages_path <- file.path(project_path, "data", "ewe", "7ages")
sa_path <- file.path(project_path, "data", "AtlanticMenhadenSA")
file_names <- c(
"biomass_annual.csv",
"catch_annual.csv",
"consumption-biomass_annual.csv",
"feedingtime_annual.csv",
"mortality_annual.csv",
"tl_annual.csv",
"weight_annual.csv"
)
name_vec <- c("biomass", "catch", "consumptionBiomass", "feedingTime", "mortality", "tl", "weight")
ewe2ages <- IP4EBFM::read_ewe_output(
file_path = file.path(ewe2ages_path, "ecosim_output"),
file_names = file_names,
skip_nrows = 8,
plot = FALSE,
figure_titles = NULL,
functional_groups = c(
"StripedBass0",
"StripedBass2_5",
"StripedBass6",
"0",
"total_adult",
"SpinyDogfish",
"BluefishJuvenile",
"BluefishAdult",
"WeakfishJuvenile",
"WeakfishAdult",
"AtlanticHerring0_1",
"AtlanticHerring2",
"Anchovies",
"Benthos",
"Zooplankton",
"Phytoplankton",
"Detritus"
),
figure_colors = NULL
)
names(ewe2ages) <- name_vec
ewe2agesAM <- c("X0", "total_adult")
ewe7agesAM <- c(
"X0",
"X1",
"X2",
"X3",
"X4",
"X5",
"X6"
)
ewe7ages_BA0 <- IP4EBFM::read_ewe_output(
file_path = file.path(ewe7ages_path, "ecosim_output", "BA0"),
file_names = file_names,
skip_nrows = 8,
plot = FALSE,
figure_titles = NULL,
functional_groups = c(
"StripedBass0",
"StripedBass2_5",
"StripedBass6",
"0",
"1",
"2",
"3",
"4",
"5",
"6",
"SpinyDogfish",
"BluefishJuvenile",
"BluefishAdult",
"WeakfishJuvenile",
"WeakfishAdult",
"AtlanticHerring0_1",
"AtlanticHerring2",
"Anchovies",
"Benthos",
"Zooplankton",
"Phytoplankton",
"Detritus"
),
figure_colors = NULL
)
names(ewe7ages_BA0) <- name_vec
ewe7ages_BAneg <- IP4EBFM::read_ewe_output(
file_path = file.path(ewe7ages_path, "ecosim_output", "BAneg"),
file_names = file_names,
skip_nrows = 8,
plot = FALSE,
figure_titles = NULL,
functional_groups = c(
"StripedBass0",
"StripedBass2_5",
"StripedBass6",
"0",
"1",
"2",
"3",
"4",
"5",
"6",
"SpinyDogfish",
"BluefishJuvenile",
"BluefishAdult",
"WeakfishJuvenile",
"WeakfishAdult",
"AtlanticHerring0_1",
"AtlanticHerring2",
"Anchovies",
"Benthos",
"Zooplankton",
"Phytoplankton",
"Detritus"
),
figure_colors = NULL
)
names(ewe7ages_BAneg) <- name_vec
ewe7ages_BApos <- IP4EBFM::read_ewe_output(
file_path = file.path(ewe7ages_path, "ecosim_output", "BApos"),
file_names = file_names,
skip_nrows = 8,
plot = FALSE,
figure_titles = NULL,
functional_groups = c(
"StripedBass0",
"StripedBass2_5",
"StripedBass6",
"0",
"1",
"2",
"3",
"4",
"5",
"6",
"SpinyDogfish",
"BluefishJuvenile",
"BluefishAdult",
"WeakfishJuvenile",
"WeakfishAdult",
"AtlanticHerring0_1",
"AtlanticHerring2",
"Anchovies",
"Benthos",
"Zooplankton",
"Phytoplankton",
"Detritus"
),
figure_colors = NULL
)
names(ewe7ages_BApos) <- name_vec
sa_data <- dget(file.path(sa_path, "am019.rdat"))
# Create a list of data
data_name <- c("SA", "EwE AM 2 ages", "EwE AM 7 ages_BA0", "EwE AM 7 ages_BAneg")
data <- vector(mode = "list", length = length(data_name))
names(data) <- data_name
sim_years <- 1985:2017
full_ages <- 0:6
adult_ages <- 1:6
sa_name_vec <- paste0("X", full_ages)
data$SA$baa <- as.data.frame(sa_data$B.mdyr.age[as.character(sim_years), as.character(full_ages)] / 1000)
names(data$SA$baa) <- sa_name_vec
data$SA$baa$total_adult <- apply(data$SA$baa[, sa_name_vec], 1, sum)
data$SA$zaa <- as.data.frame((sa_data$Z.age[as.character(sim_years), as.character(full_ages)]))
names(data$SA$zaa) <- sa_name_vec
data$SA$zaa$total_adult <- apply(data$SA$zaa[, sa_name_vec], 1, max) # max Z
data$SA$caa <- as.data.frame(sa_data$L.age.pred.mt[as.character(sim_years), as.character(full_ages)] / 1000)
names(data$SA$caa) <- sa_name_vec
data$SA$caa$total_adult <- apply(data$SA$caa[, sa_name_vec], 1, sum)
data$`EwE AM 2 ages`$baa <- ewe2ages$biomass[, ewe2agesAM]
data$`EwE AM 2 ages`$zaa <- ewe2ages$mortality[, ewe2agesAM]
data$`EwE AM 2 ages`$caa <- ewe2ages$catch[, ewe2agesAM]
data$`EwE AM 7 ages_BA0`$baa <- ewe7ages_BA0$biomass[, ewe7agesAM]
data$`EwE AM 7 ages_BA0`$baa$total_adult <- apply(data$`EwE AM 7 ages_BA0`$baa[, ewe7agesAM], 1, sum)
data$`EwE AM 7 ages_BA0`$zaa <- ewe7ages_BA0$mortality[, ewe7agesAM]
data$`EwE AM 7 ages_BA0`$zaa$total_adult <- apply(data$`EwE AM 7 ages_BA0`$zaa[, ewe7agesAM], 1, max)
data$`EwE AM 7 ages_BA0`$caa <- ewe7ages_BA0$catch[, ewe7agesAM]
data$`EwE AM 7 ages_BA0`$caa$total_adult <- apply(data$`EwE AM 7 ages_BA0`$caa[, ewe7agesAM], 1, sum)
data$`EwE AM 7 ages_BAneg`$baa <- ewe7ages_BAneg$biomass[, ewe7agesAM]
data$`EwE AM 7 ages_BAneg`$baa$total_adult <- apply(data$`EwE AM 7 ages_BAneg`$baa[, ewe7agesAM], 1, sum)
data$`EwE AM 7 ages_BAneg`$zaa <- ewe7ages_BAneg$mortality[, ewe7agesAM]
data$`EwE AM 7 ages_BAneg`$zaa$total_adult <- apply(data$`EwE AM 7 ages_BAneg`$zaa[, ewe7agesAM], 1, max)
data$`EwE AM 7 ages_BAneg`$caa <- ewe7ages_BAneg$catch[, ewe7agesAM]
data$`EwE AM 7 ages_BAneg`$caa$total_adult <- apply(data$`EwE AM 7 ages_BAneg`$caa[, ewe7agesAM], 1, sum)
# Make comparison figures -------------------------------------------------
pdf(file = file.path(project_path, "data", "ewe", "bam_ewe_comparison.pdf"), onefile = T)
par(mfrow = c(3, 3))
age_labels <- paste("Age", c("0", "1", "2", "3", "4", "5", "6+", "1-6+"))
col_names <- c(sa_name_vec, "total_adult")
# Biomass
for (i in 1:length(age_labels)) {
plot(sim_years, data$SA$baa[, i], xlab = "Year", ylab = paste("Biomass", age_labels[i]), ylim = c(0, max(data$SA$baa[, i])))
if (i %in% c(1, length(age_labels))) lines(sim_years, data$`EwE AM 2 ages`$baa[, col_names[i]], pch = 2, col = 2, type = "o")
lines(sim_years, data$`EwE AM 7 ages_BA0`$baa[, col_names[i]], pch = 3, col = 3, type = "o")
lines(sim_years, data$`EwE AM 7 ages_BAneg`$baa[, col_names[i]], pch = 4, col = 4, type = "o")
}
plot.new()
legend("bottomright",
data_name,
col = 1:4,
pch = 1:4,
lty = c(NA, 1, 1, 1),
bty = "n"
)
# Mortality
for (i in 1:length(age_labels)) {
plot(sim_years, data$SA$zaa[, i], xlab = "Year", ylab = paste("Mortality", age_labels[i]), ylim = c(0, 2))
if (i %in% c(1, length(age_labels))) lines(sim_years, data$`EwE AM 2 ages`$zaa[, col_names[i]], pch = 2, col = 2, type = "o")
lines(sim_years, data$`EwE AM 7 ages_BA0`$zaa[, col_names[i]], pch = 3, col = 3, type = "o")
lines(sim_years, data$`EwE AM 7 ages_BAneg`$zaa[, col_names[i]], pch = 4, col = 4, type = "o")
}
plot.new()
legend("bottomright",
data_name,
col = 1:4,
pch = 1:4,
lty = c(NA, 1, 1, 1),
bty = "n"
)
# Catch
for (i in 1:length(age_labels)) {
plot(sim_years, data$SA$caa[, i], xlab = "Year", ylab = paste("Catch", age_labels[i]), ylim = c(0, max(data$SA$caa[, i])))
if (i %in% c(1, length(age_labels))) lines(sim_years, data$`EwE AM 2 ages`$caa[, col_names[i]], pch = 2, col = 2, type = "o")
lines(sim_years, data$`EwE AM 7 ages_BA0`$caa[, col_names[i]], pch = 3, col = 3, type = "o")
lines(sim_years, data$`EwE AM 7 ages_BAneg`$caa[, col_names[i]], pch = 4, col = 4, type = "o")
}
plot.new()
legend("bottomright",
data_name,
col = 1:4,
pch = 1:4,
lty = c(NA, 1, 1, 1),
bty = "n"
)
dev.off()
# Biomass and catch trends from EwE of other species
library(tidyverse)
library(reshape2)
labels_2ages <- c(
"Striped Bass 0", "Striped Bass 2-5", "Striped Bass 6+",
"Atlantic Menhaden 0", "Atlantic Menhaden 1-6+",
"Spiny Dogfish", "Bluefish juv", "Bluefish adult",
"Weakfish juv", "Weakfish adult",
"Atlantic Herring 0-1", "Atlantic Herring 2+",
"Anchovies", "Benthos", "Zooplkanton", "Phytoplankton", "Detritus"
)
labels_7ages <- c(
"Striped Bass 0", "Striped Bass 2-5", "Striped Bass 6+",
"Atlantic Menhaden 0", "Atlantic Menhaden 1",
"Atlantic Menhaden 2", "Atlantic Menhaden 3",
"Atlantic Menhaden 4", "Atlantic Menhaden 5",
"Atlantic Menhaden 6+",
"Spiny Dogfish", "Bluefish juv", "Bluefish adult",
"Weakfish juv", "Weakfish adult",
"Atlantic Herring 0-1", "Atlantic Herring 2+",
"Anchovies", "Benthos", "Zooplkanton", "Phytoplankton", "Detritus",
"Atlantic Menhaden 1-6+"
)
names(ewe2ages$biomass)[2:ncol(ewe2ages$biomass)] <- labels_2ages
ewe2ages_melt <- reshape2::melt(ewe2ages$biomass, id.vars = "Year")
ewe2ages_melt$Model <- "2 age classes"
ewe7ages <- ewe7ages_BA0
ewe7ages$biomass$total_adult <- apply(ewe7ages$biomass[, ewe7agesAM], 1, sum)
names(ewe7ages$biomass)[2:ncol(ewe7ages$biomass)] <- labels_7ages
ewe7ages_melt <- reshape2::melt(ewe7ages$biomass, id.vars = "Year")
ewe7ages_melt$Model <- "7 age classes"
ewe_data <- rbind(ewe2ages_melt, ewe7ages_melt)
ggplot(ewe_data, aes(x = Year, y = value, color = Model)) +
geom_line() +
facet_wrap(~variable, scales = "free_y") +
labs(
title = "Comparison between two-ages EwE model and seven-ages EwE model",
x = "Year",
y = "Biomass (million mt)"
) +
theme_bw()
## Catch
names(ewe2ages$catch)[2:ncol(ewe2ages$catch)] <- labels_2ages
ewe2ages_melt <- reshape2::melt(ewe2ages$catch, id.vars = "Year")
ewe2ages_melt$Model <- "2 age classes"
ewe7ages <- ewe7ages_BA0
ewe7ages$catch$total_adult <- apply(ewe7ages$catch[, ewe7agesAM], 1, sum)
names(ewe7ages$catch)[2:ncol(ewe7ages$catch)] <- labels_7ages
ewe7ages_melt <- reshape2::melt(ewe7ages$catch, id.vars = "Year")
ewe7ages_melt$Model <- "7 age classes"
ewe_data <- rbind(ewe2ages_melt, ewe7ages_melt)
ggplot(ewe_data, aes(x = Year, y = value, color = Model)) +
geom_line() +
facet_wrap(~variable, scales = "free_y") +
labs(
title = "Comparison between two-ages EwE model and seven-ages EwE model",
x = "Year",
y = "Catch (million mt)"
) +
theme_bw()
Biomass-at-age of menhaden-like species over years (1985 - 2017)
``` {r, echo=FALSE, results='hide', message=FALSE, warning=FALSE, fig.width=12, fig.height=10}
plot_ewe_agecomp <- function(ewe, sa=NULL, title, xlab, ylab) {
ewe_melt <- reshape2::melt(ewe, id.vars = "Year")
ewe_melt$variable <- as.factor(ewe_melt$variable)
colnames(ewe_melt)[which(colnames(ewe_melt) == "variable")] <- "Age"
if(!is.null(sa)) {
sa_melt <- reshape2::melt(sa, id.vars = "Year")
sa_melt$variable <- as.factor(sa_melt$variable)
colnames(sa_melt)[which(colnames(sa_melt) == "variable")] <- "Age"
ggplot(data=ewe_melt, aes(x=Age, fill=Age, y=value)) +
geom_bar(position='dodge', stat='identity') +
geom_point(data=sa_melt, mapping=aes(x=Age, y=value)) +
facet_wrap(~Year) +
labs(title = title,
x = xlab,
y = ylab) +
theme_bw()
} else{
ggplot(data=ewe_melt, aes(x=Age, fill=Age, y=value)) +
geom_bar(position='dodge', stat='identity') +
facet_wrap(~Year) +
labs(title = title,
x = xlab,
y = ylab) +
theme_bw()
}
}
ewe <- ewe7ages_BA0$biomass[, c("Year", paste0("X", 0:6))]
ewe[, paste0("X", 0:6)] <- ewe[, paste0("X", 0:6)]/rowSums(ewe[, paste0("X", 0:6)])
colnames(ewe) <- c("Year", paste(0:5), "6+")
sa <- as.data.frame(sa_data$B.mdyr.age[as.numeric(row.names(sa_data$B.mdyr.age)) %in% sim_years, ]/1000)
sa <- sa/rowSums(sa)
colnames(sa) <- c(paste(0:5), "6+")
sa$Year <- sim_years
plot_ewe_agecomp(ewe=ewe,
sa=sa,
title="Biomass-at-age",
xlab="Age",
ylab = "Proportion")
The bars represent the outputs from EwE and the solid points represent the outputs from BAM stock assessment outputs.
## Mortality-at-age of menhaden-like species over years (1985 - 2017)
``` {r, echo=FALSE, results='hide', message=FALSE, warning=FALSE, fig.width=12, fig.height=10}
ewe <- ewe7ages_BA0$mortality[, c("Year", paste0("X", 0:6))]
colnames(ewe) <- c("Year", paste(0:5), "6+")
sa <- as.data.frame(sa_data$Z.age[as.numeric(row.names(sa_data$Z.age)) %in% sim_years, ])
colnames(sa) <- c(paste(0:5), "6+")
sa$Year <- sim_years
plot_ewe_agecomp(ewe=ewe,
sa=sa,
title="Mortality-at-age",
xlab="Age",
ylab = "Mortality")
Catch-at-age of menhaden-like species over years (1985 - 2017)
``` {r, echo=FALSE, results='hide', message=FALSE, warning=FALSE, fig.width=12, fig.height=10}
ewe <- ewe7ages_BA0$catch[, c("Year", paste0("X", 0:6))]
ewe[, paste0("X", 0:6)] <- ewe[, paste0("X", 0:6)]/rowSums(ewe[, paste0("X", 0:6)])
colnames(ewe) <- c("Year", paste(0:5), "6+")
sa <- as.data.frame(sa_data$L.age.pred.mt[as.numeric(row.names(sa_data$L.age.pred.mt)) %in% sim_years, ]/1000)
sa <- sa/rowSums(sa)
colnames(sa) <- c(paste(0:5), "6+")
sa$Year <- sim_years
plot_ewe_agecomp(ewe=ewe,
sa=sa,
title="Catch-at-age",
xlab="Age",
ylab = "Proportion")
## Weight-at-age of menhaden-like spcies over years (1985 - 2017)
``` {r, echo=FALSE, results='hide', message=FALSE, warning=FALSE, fig.width=12, fig.height=10}
ewe <- ewe7ages_BA0$weight[, c("Year", paste0("X", 0:6))]
colnames(ewe) <- c("Year", paste(0:5), "6+")
sa <- as.data.frame(do.call("rbind", replicate(
length(sim_years), sa_data$a.series$weight.spawn*1000, simplify = FALSE)))
colnames(sa) <- c(paste(0:5), "6+")
sa$Year <- sim_years
plot_ewe_agecomp(ewe=ewe,
sa=sa,
title="Weight-at-age",
xlab="Age",
ylab = "Weight")
Acknowledgements
Thanks Sarah Gaichas, Howard Townsend, Desiree Tommasi, and Isaac Kaplan for helping lay out the tasks to restructure the EwE model and discussing the diagnostic process.
Thanks David Chagaris and Amy Schueller for sharing model input files and providing suggestions to restructure the EwE model.
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 = TRUE)
The Ecopath with Ecosim (EwE) Model
The goal of the this work is to restructure an EwE model that include a menhaden-like species with 7 age classes (i.e. 0 - 6+). An existing EwE model of the Northwest Atlantic Continental Shelf model of intermediate complexity for ecosystems (NWACS-MICE) was previously developed to link the dynamics of menhaden with key managed predators (Chagaris et al., 2020). There were two age classes of menhaden in the NWACS-MICE model. To restructure the MWACS-MICE model, we made following changes:
-
Updated Ecopath multi-stanza groups (menhaden-like species)
-
Updated age 6+ biomass in 1985 using Beaufort Assessment Model (BAM) outputs (SEDAR 2020).
-
Updated total mortality-at-age in 1985 inputs using $C_{a}/B_{a}+M_{a}$ with BAM outputs. Here $C$, $B$, and $M$ represent catch, biomass, and natural mortality respectively. $a$ denotes age.
-
Updated relative biomass accumulation rate (BA/B). Three scenarios were explored: 1) 0.114 from Chagaris et al. (2020), 2) EwE default value of 0, 3) -0.081 based on $(B_{1986} - B_{1985})/B_{1986}$ (Buchheister et al., 2017). Both BA/B = 0 and BA/B = -0.081 help improve agreement between Ecosim predicted catch with observed catch and BA/B = 0 showed lower sum of squares (SS) value. SS values from scenarios 1 - 3 were 3059, 2339, and 3136 respectively. So we used BA/B = 0 for base run. Key figures from the three scenarios could be found here.
-
Updated Ecopath diet composition using inputs from NWACS-MICE model (Chagaris et al., 2020), full NWACS model (Buchhester et al., 2017), and BAM M-at-age data.
-
Updated Ecopath C-at-age using BAM outputs.
-
Updated Ecosim vulnerabilities and ages 1-6+ share the same vulnerability values.
-
Updated time series matrix using I-at-age, C-at-age, and F-at-age from BAM outputs. Here $I$ and $F$ represent abundance indices and fishing mortality respectively.
-
A few forcing scenarios were tested to explore various sources of environmental variability. For example, primary production splines was applied to phytoplankton, recruitment deviations from Chagaris et al. (2020) and BAM outputs were linked with menhaden age 0 fish. However, incorporating primary production splines did help soak up some of the process variability in the food web dynamics and did not help improve model fits. Incorporation of recruitment deviations did not help better represent interannual variability in the system. Question: Maybe not apply forcing functions to the EwE model, but exploring relationships between environmental drivers and menhaden biomass when adjusting $F_{ECO}$?.
-
Density-dependent growth: According to EwE 6.0 guide, users could set a non-zero feeding time adjustment for the juvenile group and this represents density-dependent changes in juvenile mortality rate associated with changes in feeding time and predation risk. The feeding time adjust rate for age 0 of menhaden-like species was set to 0.5.
Comparison between EwE model with 2 age classes of menhaden (Chagaris et al., 2020) and EwE model with 7 age classes of menhaden-like species.
devtools::load_all() # Load output data from BAM and EwE --------------------------------------- project_path <- here::here() ewe2ages_path <- file.path(project_path, "data", "ewe", "2ages") ewe7ages_path <- file.path(project_path, "data", "ewe", "7ages") sa_path <- file.path(project_path, "data", "AtlanticMenhadenSA") file_names <- c( "biomass_annual.csv", "catch_annual.csv", "consumption-biomass_annual.csv", "feedingtime_annual.csv", "mortality_annual.csv", "tl_annual.csv", "weight_annual.csv" ) name_vec <- c("biomass", "catch", "consumptionBiomass", "feedingTime", "mortality", "tl", "weight") ewe2ages <- IP4EBFM::read_ewe_output( file_path = file.path(ewe2ages_path, "ecosim_output"), file_names = file_names, skip_nrows = 8, plot = FALSE, figure_titles = NULL, functional_groups = c( "StripedBass0", "StripedBass2_5", "StripedBass6", "0", "total_adult", "SpinyDogfish", "BluefishJuvenile", "BluefishAdult", "WeakfishJuvenile", "WeakfishAdult", "AtlanticHerring0_1", "AtlanticHerring2", "Anchovies", "Benthos", "Zooplankton", "Phytoplankton", "Detritus" ), figure_colors = NULL ) names(ewe2ages) <- name_vec ewe2agesAM <- c("X0", "total_adult") ewe7agesAM <- c( "X0", "X1", "X2", "X3", "X4", "X5", "X6" ) ewe7ages_BA0 <- IP4EBFM::read_ewe_output( file_path = file.path(ewe7ages_path, "ecosim_output", "BA0"), file_names = file_names, skip_nrows = 8, plot = FALSE, figure_titles = NULL, functional_groups = c( "StripedBass0", "StripedBass2_5", "StripedBass6", "0", "1", "2", "3", "4", "5", "6", "SpinyDogfish", "BluefishJuvenile", "BluefishAdult", "WeakfishJuvenile", "WeakfishAdult", "AtlanticHerring0_1", "AtlanticHerring2", "Anchovies", "Benthos", "Zooplankton", "Phytoplankton", "Detritus" ), figure_colors = NULL ) names(ewe7ages_BA0) <- name_vec ewe7ages_BAneg <- IP4EBFM::read_ewe_output( file_path = file.path(ewe7ages_path, "ecosim_output", "BAneg"), file_names = file_names, skip_nrows = 8, plot = FALSE, figure_titles = NULL, functional_groups = c( "StripedBass0", "StripedBass2_5", "StripedBass6", "0", "1", "2", "3", "4", "5", "6", "SpinyDogfish", "BluefishJuvenile", "BluefishAdult", "WeakfishJuvenile", "WeakfishAdult", "AtlanticHerring0_1", "AtlanticHerring2", "Anchovies", "Benthos", "Zooplankton", "Phytoplankton", "Detritus" ), figure_colors = NULL ) names(ewe7ages_BAneg) <- name_vec ewe7ages_BApos <- IP4EBFM::read_ewe_output( file_path = file.path(ewe7ages_path, "ecosim_output", "BApos"), file_names = file_names, skip_nrows = 8, plot = FALSE, figure_titles = NULL, functional_groups = c( "StripedBass0", "StripedBass2_5", "StripedBass6", "0", "1", "2", "3", "4", "5", "6", "SpinyDogfish", "BluefishJuvenile", "BluefishAdult", "WeakfishJuvenile", "WeakfishAdult", "AtlanticHerring0_1", "AtlanticHerring2", "Anchovies", "Benthos", "Zooplankton", "Phytoplankton", "Detritus" ), figure_colors = NULL ) names(ewe7ages_BApos) <- name_vec sa_data <- dget(file.path(sa_path, "am019.rdat")) # Create a list of data data_name <- c("SA", "EwE AM 2 ages", "EwE AM 7 ages_BA0", "EwE AM 7 ages_BAneg") data <- vector(mode = "list", length = length(data_name)) names(data) <- data_name sim_years <- 1985:2017 full_ages <- 0:6 adult_ages <- 1:6 sa_name_vec <- paste0("X", full_ages) data$SA$baa <- as.data.frame(sa_data$B.mdyr.age[as.character(sim_years), as.character(full_ages)] / 1000) names(data$SA$baa) <- sa_name_vec data$SA$baa$total_adult <- apply(data$SA$baa[, sa_name_vec], 1, sum) data$SA$zaa <- as.data.frame((sa_data$Z.age[as.character(sim_years), as.character(full_ages)])) names(data$SA$zaa) <- sa_name_vec data$SA$zaa$total_adult <- apply(data$SA$zaa[, sa_name_vec], 1, max) # max Z data$SA$caa <- as.data.frame(sa_data$L.age.pred.mt[as.character(sim_years), as.character(full_ages)] / 1000) names(data$SA$caa) <- sa_name_vec data$SA$caa$total_adult <- apply(data$SA$caa[, sa_name_vec], 1, sum) data$`EwE AM 2 ages`$baa <- ewe2ages$biomass[, ewe2agesAM] data$`EwE AM 2 ages`$zaa <- ewe2ages$mortality[, ewe2agesAM] data$`EwE AM 2 ages`$caa <- ewe2ages$catch[, ewe2agesAM] data$`EwE AM 7 ages_BA0`$baa <- ewe7ages_BA0$biomass[, ewe7agesAM] data$`EwE AM 7 ages_BA0`$baa$total_adult <- apply(data$`EwE AM 7 ages_BA0`$baa[, ewe7agesAM], 1, sum) data$`EwE AM 7 ages_BA0`$zaa <- ewe7ages_BA0$mortality[, ewe7agesAM] data$`EwE AM 7 ages_BA0`$zaa$total_adult <- apply(data$`EwE AM 7 ages_BA0`$zaa[, ewe7agesAM], 1, max) data$`EwE AM 7 ages_BA0`$caa <- ewe7ages_BA0$catch[, ewe7agesAM] data$`EwE AM 7 ages_BA0`$caa$total_adult <- apply(data$`EwE AM 7 ages_BA0`$caa[, ewe7agesAM], 1, sum) data$`EwE AM 7 ages_BAneg`$baa <- ewe7ages_BAneg$biomass[, ewe7agesAM] data$`EwE AM 7 ages_BAneg`$baa$total_adult <- apply(data$`EwE AM 7 ages_BAneg`$baa[, ewe7agesAM], 1, sum) data$`EwE AM 7 ages_BAneg`$zaa <- ewe7ages_BAneg$mortality[, ewe7agesAM] data$`EwE AM 7 ages_BAneg`$zaa$total_adult <- apply(data$`EwE AM 7 ages_BAneg`$zaa[, ewe7agesAM], 1, max) data$`EwE AM 7 ages_BAneg`$caa <- ewe7ages_BAneg$catch[, ewe7agesAM] data$`EwE AM 7 ages_BAneg`$caa$total_adult <- apply(data$`EwE AM 7 ages_BAneg`$caa[, ewe7agesAM], 1, sum) # Make comparison figures ------------------------------------------------- pdf(file = file.path(project_path, "data", "ewe", "bam_ewe_comparison.pdf"), onefile = T) par(mfrow = c(3, 3)) age_labels <- paste("Age", c("0", "1", "2", "3", "4", "5", "6+", "1-6+")) col_names <- c(sa_name_vec, "total_adult") # Biomass for (i in 1:length(age_labels)) { plot(sim_years, data$SA$baa[, i], xlab = "Year", ylab = paste("Biomass", age_labels[i]), ylim = c(0, max(data$SA$baa[, i]))) if (i %in% c(1, length(age_labels))) lines(sim_years, data$`EwE AM 2 ages`$baa[, col_names[i]], pch = 2, col = 2, type = "o") lines(sim_years, data$`EwE AM 7 ages_BA0`$baa[, col_names[i]], pch = 3, col = 3, type = "o") lines(sim_years, data$`EwE AM 7 ages_BAneg`$baa[, col_names[i]], pch = 4, col = 4, type = "o") } plot.new() legend("bottomright", data_name, col = 1:4, pch = 1:4, lty = c(NA, 1, 1, 1), bty = "n" ) # Mortality for (i in 1:length(age_labels)) { plot(sim_years, data$SA$zaa[, i], xlab = "Year", ylab = paste("Mortality", age_labels[i]), ylim = c(0, 2)) if (i %in% c(1, length(age_labels))) lines(sim_years, data$`EwE AM 2 ages`$zaa[, col_names[i]], pch = 2, col = 2, type = "o") lines(sim_years, data$`EwE AM 7 ages_BA0`$zaa[, col_names[i]], pch = 3, col = 3, type = "o") lines(sim_years, data$`EwE AM 7 ages_BAneg`$zaa[, col_names[i]], pch = 4, col = 4, type = "o") } plot.new() legend("bottomright", data_name, col = 1:4, pch = 1:4, lty = c(NA, 1, 1, 1), bty = "n" ) # Catch for (i in 1:length(age_labels)) { plot(sim_years, data$SA$caa[, i], xlab = "Year", ylab = paste("Catch", age_labels[i]), ylim = c(0, max(data$SA$caa[, i]))) if (i %in% c(1, length(age_labels))) lines(sim_years, data$`EwE AM 2 ages`$caa[, col_names[i]], pch = 2, col = 2, type = "o") lines(sim_years, data$`EwE AM 7 ages_BA0`$caa[, col_names[i]], pch = 3, col = 3, type = "o") lines(sim_years, data$`EwE AM 7 ages_BAneg`$caa[, col_names[i]], pch = 4, col = 4, type = "o") } plot.new() legend("bottomright", data_name, col = 1:4, pch = 1:4, lty = c(NA, 1, 1, 1), bty = "n" ) dev.off() # Biomass and catch trends from EwE of other species library(tidyverse) library(reshape2) labels_2ages <- c( "Striped Bass 0", "Striped Bass 2-5", "Striped Bass 6+", "Atlantic Menhaden 0", "Atlantic Menhaden 1-6+", "Spiny Dogfish", "Bluefish juv", "Bluefish adult", "Weakfish juv", "Weakfish adult", "Atlantic Herring 0-1", "Atlantic Herring 2+", "Anchovies", "Benthos", "Zooplkanton", "Phytoplankton", "Detritus" ) labels_7ages <- c( "Striped Bass 0", "Striped Bass 2-5", "Striped Bass 6+", "Atlantic Menhaden 0", "Atlantic Menhaden 1", "Atlantic Menhaden 2", "Atlantic Menhaden 3", "Atlantic Menhaden 4", "Atlantic Menhaden 5", "Atlantic Menhaden 6+", "Spiny Dogfish", "Bluefish juv", "Bluefish adult", "Weakfish juv", "Weakfish adult", "Atlantic Herring 0-1", "Atlantic Herring 2+", "Anchovies", "Benthos", "Zooplkanton", "Phytoplankton", "Detritus", "Atlantic Menhaden 1-6+" ) names(ewe2ages$biomass)[2:ncol(ewe2ages$biomass)] <- labels_2ages ewe2ages_melt <- reshape2::melt(ewe2ages$biomass, id.vars = "Year") ewe2ages_melt$Model <- "2 age classes" ewe7ages <- ewe7ages_BA0 ewe7ages$biomass$total_adult <- apply(ewe7ages$biomass[, ewe7agesAM], 1, sum) names(ewe7ages$biomass)[2:ncol(ewe7ages$biomass)] <- labels_7ages ewe7ages_melt <- reshape2::melt(ewe7ages$biomass, id.vars = "Year") ewe7ages_melt$Model <- "7 age classes" ewe_data <- rbind(ewe2ages_melt, ewe7ages_melt) ggplot(ewe_data, aes(x = Year, y = value, color = Model)) + geom_line() + facet_wrap(~variable, scales = "free_y") + labs( title = "Comparison between two-ages EwE model and seven-ages EwE model", x = "Year", y = "Biomass (million mt)" ) + theme_bw() ## Catch names(ewe2ages$catch)[2:ncol(ewe2ages$catch)] <- labels_2ages ewe2ages_melt <- reshape2::melt(ewe2ages$catch, id.vars = "Year") ewe2ages_melt$Model <- "2 age classes" ewe7ages <- ewe7ages_BA0 ewe7ages$catch$total_adult <- apply(ewe7ages$catch[, ewe7agesAM], 1, sum) names(ewe7ages$catch)[2:ncol(ewe7ages$catch)] <- labels_7ages ewe7ages_melt <- reshape2::melt(ewe7ages$catch, id.vars = "Year") ewe7ages_melt$Model <- "7 age classes" ewe_data <- rbind(ewe2ages_melt, ewe7ages_melt) ggplot(ewe_data, aes(x = Year, y = value, color = Model)) + geom_line() + facet_wrap(~variable, scales = "free_y") + labs( title = "Comparison between two-ages EwE model and seven-ages EwE model", x = "Year", y = "Catch (million mt)" ) + theme_bw()
Biomass-at-age of menhaden-like species over years (1985 - 2017)
``` {r, echo=FALSE, results='hide', message=FALSE, warning=FALSE, fig.width=12, fig.height=10} plot_ewe_agecomp <- function(ewe, sa=NULL, title, xlab, ylab) { ewe_melt <- reshape2::melt(ewe, id.vars = "Year") ewe_melt$variable <- as.factor(ewe_melt$variable) colnames(ewe_melt)[which(colnames(ewe_melt) == "variable")] <- "Age"
if(!is.null(sa)) { sa_melt <- reshape2::melt(sa, id.vars = "Year") sa_melt$variable <- as.factor(sa_melt$variable) colnames(sa_melt)[which(colnames(sa_melt) == "variable")] <- "Age"
ggplot(data=ewe_melt, aes(x=Age, fill=Age, y=value)) +
geom_bar(position='dodge', stat='identity') +
geom_point(data=sa_melt, mapping=aes(x=Age, y=value)) +
facet_wrap(~Year) +
labs(title = title,
x = xlab,
y = ylab) +
theme_bw()
} else{ ggplot(data=ewe_melt, aes(x=Age, fill=Age, y=value)) + geom_bar(position='dodge', stat='identity') + facet_wrap(~Year) + labs(title = title, x = xlab, y = ylab) + theme_bw() }
}
ewe <- ewe7ages_BA0$biomass[, c("Year", paste0("X", 0:6))] ewe[, paste0("X", 0:6)] <- ewe[, paste0("X", 0:6)]/rowSums(ewe[, paste0("X", 0:6)]) colnames(ewe) <- c("Year", paste(0:5), "6+")
sa <- as.data.frame(sa_data$B.mdyr.age[as.numeric(row.names(sa_data$B.mdyr.age)) %in% sim_years, ]/1000) sa <- sa/rowSums(sa) colnames(sa) <- c(paste(0:5), "6+") sa$Year <- sim_years
plot_ewe_agecomp(ewe=ewe, sa=sa, title="Biomass-at-age", xlab="Age", ylab = "Proportion")
The bars represent the outputs from EwE and the solid points represent the outputs from BAM stock assessment outputs. ## Mortality-at-age of menhaden-like species over years (1985 - 2017) ``` {r, echo=FALSE, results='hide', message=FALSE, warning=FALSE, fig.width=12, fig.height=10} ewe <- ewe7ages_BA0$mortality[, c("Year", paste0("X", 0:6))] colnames(ewe) <- c("Year", paste(0:5), "6+") sa <- as.data.frame(sa_data$Z.age[as.numeric(row.names(sa_data$Z.age)) %in% sim_years, ]) colnames(sa) <- c(paste(0:5), "6+") sa$Year <- sim_years plot_ewe_agecomp(ewe=ewe, sa=sa, title="Mortality-at-age", xlab="Age", ylab = "Mortality")
Catch-at-age of menhaden-like species over years (1985 - 2017)
``` {r, echo=FALSE, results='hide', message=FALSE, warning=FALSE, fig.width=12, fig.height=10}
ewe <- ewe7ages_BA0$catch[, c("Year", paste0("X", 0:6))] ewe[, paste0("X", 0:6)] <- ewe[, paste0("X", 0:6)]/rowSums(ewe[, paste0("X", 0:6)]) colnames(ewe) <- c("Year", paste(0:5), "6+")
sa <- as.data.frame(sa_data$L.age.pred.mt[as.numeric(row.names(sa_data$L.age.pred.mt)) %in% sim_years, ]/1000) sa <- sa/rowSums(sa) colnames(sa) <- c(paste(0:5), "6+") sa$Year <- sim_years
plot_ewe_agecomp(ewe=ewe, sa=sa, title="Catch-at-age", xlab="Age", ylab = "Proportion")
## Weight-at-age of menhaden-like spcies over years (1985 - 2017) ``` {r, echo=FALSE, results='hide', message=FALSE, warning=FALSE, fig.width=12, fig.height=10} ewe <- ewe7ages_BA0$weight[, c("Year", paste0("X", 0:6))] colnames(ewe) <- c("Year", paste(0:5), "6+") sa <- as.data.frame(do.call("rbind", replicate( length(sim_years), sa_data$a.series$weight.spawn*1000, simplify = FALSE))) colnames(sa) <- c(paste(0:5), "6+") sa$Year <- sim_years plot_ewe_agecomp(ewe=ewe, sa=sa, title="Weight-at-age", xlab="Age", ylab = "Weight")
Acknowledgements
Thanks Sarah Gaichas, Howard Townsend, Desiree Tommasi, and Isaac Kaplan for helping lay out the tasks to restructure the EwE model and discussing the diagnostic process.
Thanks David Chagaris and Amy Schueller for sharing model input files and providing suggestions to restructure the EwE model.
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.