Multispecies Recruitment
Sarah Gaichas
08 January, 2024
# hacks for data not yet in stockSmart as of January 2024
# mackerel MT output from draft via Dan Hennen and Kiersten Curti
# see local folder mackerelMT2023 for inputs
# set addmack to FALSE when stockSmart is updated
addmack <- TRUE
# add mackerel MT results by hand
if(addmack){
mackSS <- stocksmart::stockAssessmentData |>
dplyr::filter(StockName == "Atlantic mackerel - Gulf of Maine / Cape Hatteras",
AssessmentYear == max(AssessmentYear)) |>
dplyr::select(StockName, Metric, Description, Units, AssessmentYear, Jurisdiction) |>
dplyr::group_by(Metric) |>
dplyr::distinct()
# 2013-2022 match for values in stocksmart to add units, etc from there
newmack <- read.csv(here::here("mackerelMT2023/model_results.csv")) |>
dplyr::select(Year, Abundance = SSB, Recruitment = Rect) |>
tidyr::pivot_longer(-Year, names_to = "Metric", values_to = "Value") |>
dplyr::mutate(StockName = "Atlantic mackerel - Gulf of Maine / Cape Hatteras") |>
dplyr::left_join(mackSS)
}Evaluate change in recruitment over time for assessed species
This is an update to analyses presented in Perretti et al. (2017). Code courtesy Charles Perretti and Kim Bastille.
Get stocksmart info
UPDATE: StockSMART was updated in mid-December 2023 with Northeast US
assessment results, and stocksmart now has all results.
Except! Atlantic mackerel had only 2013-2022 in stockSMART so I
am adding the draft results from the MT assessment by Kiersten Curti via
Dan Hennen.
Can we pull most recent recruitment and SSB or other B outputs from Northeast US stock assessments? Which species do we have? What units are recruitments in?
These are all stocks, most recent assessment year reporting recruitment.
NErec <- stocksmart::stockAssessmentData %>%
filter(RegionalEcosystem == "Northeast Shelf",
Metric == "Recruitment") %>%
group_by(StockName) %>%
filter(AssessmentYear == max(AssessmentYear)) %>%
ungroup()
if(addmack) {
macrec <- newmack |>
dplyr::filter(Metric == "Recruitment")
NErec <- full_join(NErec, macrec)
}
sppunits <- NErec %>%
select(StockName, StockArea, Description, Units) %>%
distinct()
datatable(sppunits, rownames = FALSE,
options = list(scrollX = TRUE,
pageLength = 10))These are ball stocks, most recent recruitment year, reporting abundance. Lets check the units for biomass too.
NEbio <- stocksmart::stockAssessmentData %>%
filter(RegionalEcosystem == "Northeast Shelf",
Metric == "Abundance") %>%
group_by(StockName) %>%
filter(AssessmentYear == max(AssessmentYear)) %>%
ungroup()
if(addmack) {
macbio <- newmack |>
dplyr::filter(Metric == "Abundance")
NEbio <- full_join(NEbio, macbio)
}
sppunits <- NEbio %>%
select(StockName, StockArea, Description, Units) %>%
distinct()
datatable(sppunits, rownames = FALSE,
options = list(scrollX = TRUE,
pageLength = 10))Select down to stocks used in SOE productivity anomaly indicator? Won’t split them by EPU but could map some stocks to different regions kinda. Or just pick MA managed stocks in stocksmart? Here just trying to get a list of stocks that have both recruitment and biomass in stocksmart. Weeding out recruitment measures that are not in numbers for now.
Also, lets not use assessments prior to 2019.
ADD A CHECK: do we only have one assessment per species? Yes, we did pull max(AssessmentYear)
Which assessments are in stocksmart with assessment 2019
or later? Spawning stock biomass on the left and recruitment on the
right:
These don’t align well. How about a useless sentence?
The input assessment set
Which stocks with assessments 2019 or newer have both recruitment in numbers and SSB in biomass? This limits us to fewer stocks.
# first pass only stocks with both rec in N and bio but keep all years of bio
NErecN <- NErec %>%
filter(AssessmentYear>2018,
Units %in% c("Thousand Recruits",
"Number x 1,000,000",
"Number x 1,000",
"Million Recruits"))
NErecstocks <- NErecN %>%
select(StockName) %>%
distinct()
NEbiohasrec <- NEbio %>%
filter(StockName %in% NErecstocks$StockName)
bothyr <- NEbiohasrec %>%
group_by(StockName) %>%
select(StockName, AssessmentYear) %>%
distinct()
datatable(bothyr)Standardize assessment outputs
Next steps:
Reconcile the units so inputs are the same across stocks
Units of recruitment in numbers of fish at age 1 is the input to the survey based productivity anomaly indicator. We’ll need to convert all to numbers and note the age of recruitment for each stock to shift recruitment appropriately.
NErecNstd <- NErecN %>%
mutate(NfishRec = case_when(Units == "Thousand Recruits" ~ Value*1000,
Units == "Number x 1,000,000" ~ Value*1000000,
Units == "Number x 1,000" ~ Value*1000,
Units == "Million Recruits" ~ Value*1000000,
TRUE ~ as.numeric(NA)),
AgeRec = as.numeric(stringr::str_extract(Description, "(\\d)+")),
SSByr = Year-AgeRec
)Plot each Rec series by SSByr
ggplot(NErecNstd, aes(x=SSByr, y=NfishRec)) +
geom_line() +
theme_bw() +
facet_wrap(~StockName, ncol = 3, scales = "free",
labeller = labeller(StockName = label_wrap_gen(35)))
Units of biomass (all listed as mature/spawning stock or retro adjusted spawning stock) are converted to kg for similarity with survey anomaly code.
In 2023 the MT assessment for spiny dogfish now reports stock size as reproductive potential in million pups. In order to keep spiny dogfish in the mix for this analysis, I am converting this to biomass units assuming that a pup weighs 0.5 kg. We are calculating an anomaly of recruits per biomass so I think this is ok even though it is not a measure of adult biomass.
NEbioBstd <- NEbiohasrec %>%
mutate(SSBkg = case_when(Units == "Thousand Metric Tons" ~ Value*1000000,
Units == "Metric Tons" ~ Value*1000,
Units == "Million Pups" ~ Value*10000000/2,
TRUE ~ as.numeric(NA))
#AgeRec = as.numeric(stringr::str_extract(Description, "(\\d)+")),
#SSByr = Year-AgeRec
)Plot each SSB series by Year
ggplot(NEbioBstd, aes(x=Year, y=SSBkg)) +
geom_line() +
theme_bw() +
facet_wrap(~StockName, ncol = 3, scales = "free",
labeller = labeller(StockName = label_wrap_gen(35)))
Do the Perretti anomaly calculation
For the survey data, this is in trawlr package
func-data.R, function calc_rec() starting line
798.
Current SOE indicator using rs_anom output of this as
input for plotting so need to calculate equivalent of that one.
Modify code to take stocksmart inputs converted to
common units. Because the survey anomalies are calculated in numbers of
small fish per kg mature biomass we will convert assessments to those
same units.
This is trawlr::calc_z used in
calc_rec()
# ---------------------------------------------------------
#' Calculate z-score of vector
#'
#'
#' @param x vector of which to calculate z-score
#'
#' @param ... other params that can be fed to mean() or sd()
#'
#' @return vector of z-score values of x
#'
#' @examples
#'
#' x_zscore <- calc_z(x = c(1:10), na.rm = T)
#'
#' @export
calc_z <- function(x, ...) {
(x - mean(x, ...))/sd(x, ...)
}This is trawlr::calc_rec() full code.
#' Calculate spawner-recruit anomalies
#'
#' Calculates a time series of spawner-recruit anomalies
#' (and other anomalies) for all species in survdat. Survdat
#' must have biological and epu data, so \emph{load_survdat()}
#' must have been called with \emph(with_bio_epu) = TRUE.
#'
#' @param survdat dataframe; from \emph{load_survdat()}
#' called with \emph{with_bio_epu} = TRUE.
#'
#' @param min_year numeric value; minimum year to include
#' in the recruitment calculations.
#'
#' @param max_year numeric value; maximum year to include
#' in the recruitment calculations.
#'
#' @param by_epu logical; if TRUE recruitment variables are
#' calculated per ecological production unit
#' (EPU).
#'
#' @param by_taxon logical; if TRUE recruitment variables are
#' calculated per ecological production unit
#' (EPU).
#'
#'
#' @return dataframe with recruitment related variables.
#'
#' @examples
#' calc_rec(survdat = survdat,
#' min_year = 1980, max_year = 2008,
#' by_epu = TRUE, by_taxon = FALSE)
#'
#' @export
calc_rec <- function(survdat,
min_year,
max_year,
by_epu,
by_taxon,
by_season = FALSE,
len_cutoff = 0.2) {
# Calculate the mean length at age-1 for species
# with age data
df_len_at_age1 <-
survdat %>%
dplyr::filter(YEAR >= (min_year - 1),
YEAR <= (max_year + 1),
!is.na(AGE),
!is.na(LENGTH)) %>%
dplyr::group_by(COMNAME, SCINAME) %>%
dplyr::do(fit_length_at_age1(AGE = .$AGE,
LENGTH = .$LENGTH))
# Calculate the number of tows in each survey season
dat_tows <-
survdat %>%
dplyr::distinct(CRUISE6, YEAR, SEASON, STATION, STRATUM) %>%
dplyr::select(CRUISE6, YEAR, SEASON) %>%
dplyr::group_by(CRUISE6, YEAR, SEASON) %>%
dplyr::summarise(n_tows = n())
dat_spec_rec <-
survdat %>%
dplyr::left_join(dat_tows) %>%
dplyr::left_join(df_len_at_age1) %>%
dplyr::group_by(COMNAME, SCINAME) %>%
dplyr::filter(EPU %in% c("MAB", "GB", "GOM")) %>%
dplyr::filter(YEAR >= (min_year - 1),
YEAR <= (max_year + 1)) %>%
#dplyr::do(trawlr::calc_lw(dat = .)) %>%
dplyr::group_by(., COMNAME) %>%
{ if (by_epu) dplyr::group_by(., EPU, add = T) else . } %>%
{ if (by_taxon) dplyr::group_by(., order, add = T) else . } %>%
{ if (by_season) dplyr::group_by(., SEASON, add = T) else . } %>%
dplyr::mutate(rank_LENGTH = dplyr::cume_dist(LENGTH),
NUMLEN = as.numeric(NUMLEN),
INDWT = as.numeric(INDWT),
WTperLENGTH = INDWT / LENGTH) %>%
# Remove NA's
dplyr::filter(!is.na(rank_LENGTH)) %>%
# If length at age1 exists use it, else use the length
# cutoff
dplyr::mutate(maturity = ifelse(is.na(length_at_age1),
ifelse(rank_LENGTH > len_cutoff,
"spawner",
"recruit"),
ifelse(LENGTH > length_at_age1,
"spawner",
"recruit"))
) %>%
dplyr::mutate(ind_sp_abund = ifelse(maturity == "spawner",
NUMLEN,
0),
ind_rec_abund = ifelse(maturity == "recruit",
NUMLEN,
0),
ind_sp_wl = ifelse(maturity == "spawner",
WTperLENGTH,
NaN),
ind_rec_wl = ifelse(maturity == "recruit",
WTperLENGTH,
NaN),
ind_sp_wl_res = ifelse(maturity == "spawner",
(INDWT - INDWT_ALLPRED),
NaN),
ind_rec_wl_res = ifelse(maturity == "recruit",
(INDWT - INDWT_ALLPRED),
NaN)) %>%
{ if (by_taxon == TRUE & by_epu == TRUE) {
dplyr::group_by(., YEAR, SEASON,
SCINAME, COMNAME, EPU, order)
} else if (by_epu == TRUE) {
dplyr::group_by(., YEAR, SEASON,
SCINAME, COMNAME, EPU)
} else {
dplyr::group_by(., YEAR, SEASON,
SCINAME, COMNAME)
}
} %>%
dplyr::summarise(n_tows = unique(n_tows),
spawners_abund = sum(ind_sp_abund, na.rm = T)/n_tows,
spawners_biom = sum(ind_sp_abund *
INDWT_PRED, na.rm = T)/n_tows,
recruits_abund = sum(ind_rec_abund, na.rm = T)/n_tows,
recruits_biom = sum(ind_rec_abund *
INDWT_PRED, na.rm = T)/n_tows,
spawner_wl = mean(ind_sp_wl, na.rm = T),
recruit_wl = mean(ind_rec_wl, na.rm = T),
spawner_wl_res = mean(ind_sp_wl_res,
na.rm = T),
recruit_wl_res = mean(ind_rec_wl_res,
na.rm = T)) %>%
# Calculate log(R/S), log(R), log(S)
dplyr::ungroup(.) %>%
{ if (by_taxon == TRUE & by_epu == TRUE) {
dplyr::arrange(., SEASON,
COMNAME, SCINAME, EPU, order, YEAR) %>%
dplyr::group_by(., SEASON,
COMNAME, SCINAME, EPU, order)
} else if (by_epu == TRUE) {
dplyr::arrange(., SEASON,
COMNAME, SCINAME, EPU, YEAR) %>%
dplyr::group_by(., SEASON,
COMNAME, SCINAME, EPU)
} else if (by_taxon == TRUE) {
dplyr::arrange(., SEASON,
COMNAME, SCINAME, order, YEAR) %>%
dplyr::group_by(., SEASON,
COMNAME, SCINAME, order)
} else {
dplyr::arrange(., SEASON,
COMNAME, SCINAME, YEAR) %>%
dplyr::group_by(., SEASON,
COMNAME, SCINAME)
}
} %>%
dplyr::mutate(recruits_biom_lead1 = dplyr::lead(recruits_biom,
n = 1),
recruits_abund_lead1 = dplyr::lead(recruits_abund,
n = 1),
rs = recruits_abund_lead1/
spawners_biom,
rs_abund = recruits_abund_lead1/
spawners_abund,
rs_biom = recruits_biom_lead1/
spawners_biom,
logr_abund = log(recruits_abund_lead1),
logr_biom = log(recruits_biom_lead1),
logs_abund = log(spawners_abund),
logs_biom = log(spawners_biom),
logrs = log(rs),
logrs_abund = log(rs_abund),
logrs_biom = log(rs_biom),
spawners_abund_lag0_anom =
trawlr::calc_z(spawners_abund,
na.rm = T),
spawners_biom_lag0_anom =
trawlr::calc_z(spawners_biom,
na.rm = T),
recruits_abund_lead1_anom =
trawlr::calc_z(recruits_abund_lead1,
na.rm = T),
recruits_biom_lead1_anom =
trawlr::calc_z(recruits_biom_lead1,
na.rm = T),
rs_anom = trawlr::calc_z(rs,
na.rm = T),
rs_abund_anom = trawlr::calc_z(rs_abund,
na.rm = T),
rs_biom_anom = trawlr::calc_z(rs_biom,
na.rm = T),
logr_abund_anom =
trawlr::calc_z(logr_abund),
logr_biom_anom =
trawlr::calc_z(logr_biom),
logs_abund_anom =
trawlr::calc_z(logs_abund, na.rm = T),
logs_biom_anom =
trawlr::calc_z(logs_biom, na.rm = T),
logrs_anom =
trawlr::calc_z(logrs, na.rm = T),
logrs_abund_anom =
trawlr::calc_z(logrs_abund, na.rm = T),
logrs_biom_anom =
trawlr::calc_z(logrs_biom, na.rm = T),
spawner_wl_anom =
trawlr::calc_z(spawner_wl),
recruit_wl_lead1 = dplyr::lead(recruit_wl, n = 1),
recruit_wl_lead1_anom =
trawlr::calc_z(recruit_wl_lead1),
spawner_wl_res_anom =
trawlr::calc_z(spawner_wl_res),
recruit_wl_res_lead1 = dplyr::lag(recruit_wl_res, n = 1),
recruit_wl_res_lead1_anom =
trawlr::calc_z(recruit_wl_res_lead1)) %>%
# Remove year before minimum and year after maximum
dplyr::filter(YEAR >= min_year,
YEAR <= max_year) %>%
#Average seasonal anomalies to get yearly anomaly
dplyr::group_by(., YEAR, SCINAME, COMNAME) %>%
{ if (by_epu) dplyr::group_by(., EPU, add = T) else . } %>%
{ if (by_taxon) dplyr::group_by(., order, add = T) else . } %>%
{ if (by_season) dplyr::group_by(., SEASON, add = T) else . } %>%
dplyr::summarise(
spawners_abund_lag0 = mean(spawners_abund,
na.rm = T),
spawners_abund_lag0_anom = mean(spawners_abund_lag0_anom,
na.rm = T),
spawners_biom_lag0 = mean(spawners_biom, na.rm = T),
spawners_biom_lag0_anom = mean(spawners_biom_lag0_anom,
na.rm = T),
recruits_abund = mean(recruits_abund, na.rm = T),
recruits_abund_lead1 = mean(recruits_abund_lead1,
na.rm = T),
recruits_abund_lead1_anom = mean(recruits_abund_lead1_anom,
na.rm = T),
recruits_biom_lead1 = mean(recruits_biom_lead1,
na.rm = T),
recruits_biom_lead1_anom = mean(recruits_biom_lead1_anom,
na.rm = T),
rs = mean(rs, na.rm = T),
rs_abund = mean(rs_abund, na.rm = T),
rs_biom = mean(rs_biom, na.rm = T),
rs_anom = mean(rs_anom, na.rm = T),
rs_abund_anom = mean(rs_abund_anom, na.rm = T),
rs_biom_anom = mean(rs_biom_anom, na.rm = T),
logr_abund_anom = mean(logr_abund_anom,
na.rm = T),
logr_biom_anom = mean(logr_biom_anom,
na.rm = T),
logs_abund_anom = mean(logs_abund_anom,
na.rm = T),
logs_biom_anom = mean(logs_biom_anom,
na.rm = T),
logrs_anom = mean(logrs_anom, na.rm = T),
logrs_abund_anom = mean(logrs_abund_anom,
na.rm = T),
logrs_biom_anom = mean(logrs_biom_anom,
na.rm = T),
spawner_wl_anom = mean(spawner_wl_anom,
na.rm = T),
recruit_wl_lead1_anom = mean(recruit_wl_lead1_anom,
na.rm = T),
spawner_wl_res_anom = mean(spawner_wl_res_anom,
na.rm = T),
recruit_wl_res_lead1_anom = mean(recruit_wl_res_lead1_anom,
na.rm = T))
return(dat_spec_rec)
}This is the bit I’ve pulled out to use with the stock assessment
recruitment and SSB series standardized above. Note leading and lagging
is not used as we have done this with the input assessment info by
aligning aged recruits with SSB year. We save the ecodata
submission at the end of this code chunk.
# dataset with each species std rec in SSByr and SSB
SSrecbio <- NErecNstd %>%
select(StockName, SSByr, NfishRec) %>%
left_join(NEbioBstd, by = c("StockName" = "StockName",
"SSByr"="Year")) %>%
select(StockName, SSByr, NfishRec, SSBkg)
min_year <- 1980
max_year <- 2021
dat_spec_rec <- SSrecbio %>%
group_by(StockName) %>%
dplyr::mutate(#recruits_biom_lead1 = dplyr::lead(recruits_biom,
# n = 1),
recruits_abund_lead1 = NfishRec, #dplyr::lead(recruits_abund,
# n = 1),
spawners_biom = SSBkg,
YEAR = SSByr,
rs = recruits_abund_lead1/
spawners_biom,
#rs_abund = recruits_abund_lead1/
# spawners_abund,
#rs_biom = recruits_biom_lead1/
# spawners_biom,
logr_abund = log(recruits_abund_lead1),
#logr_biom = log(recruits_biom_lead1),
#logs_abund = log(spawners_abund),
logs_biom = log(spawners_biom),
logrs = log(rs),
#logrs_abund = log(rs_abund),
#logrs_biom = log(rs_biom),
#spawners_abund_lag0_anom =
# ''calc_z(spawners_abund,
# na.rm = T),
spawners_biom_lag0_anom =
calc_z(spawners_biom,
na.rm = T),
recruits_abund_lead1_anom =
calc_z(recruits_abund_lead1,
na.rm = T),
#recruits_biom_lead1_anom =
# calc_z(recruits_biom_lead1,
# na.rm = T),
rs_anom = calc_z(rs,
na.rm = T),
#rs_abund_anom = calc_z(rs_abund,
# na.rm = T),
#rs_biom_anom = calc_z(rs_biom,
# na.rm = T),
logr_abund_anom =
calc_z(logr_abund),
#logr_biom_anom =
# calc_z(logr_biom),
#logs_abund_anom =
# calc_z(logs_abund, na.rm = T),
logs_biom_anom =
calc_z(logs_biom, na.rm = T),
logrs_anom =
calc_z(logrs, na.rm = T),
#logrs_abund_anom =
# calc_z(logrs_abund, na.rm = T),
#logrs_biom_anom =
# calc_z(logrs_biom, na.rm = T),
#spawner_wl_anom =
# calc_z(spawner_wl),
#recruit_wl_lead1 = dplyr::lead(recruit_wl, n = 1),
#recruit_wl_lead1_anom =
# calc_z(recruit_wl_lead1),
#spawner_wl_res_anom =
# calc_z(spawner_wl_res),
#recruit_wl_res_lead1 = dplyr::lag(recruit_wl_res, n = 1),
#recruit_wl_res_lead1_anom =
# calc_z(recruit_wl_res_lead1)
) %>%
# Remove year before minimum and year after maximum
dplyr::filter(YEAR >= min_year,
YEAR <= max_year) %>%
#Average seasonal anomalies to get yearly anomaly
dplyr::group_by(., YEAR, StockName) %>% #SCINAME, COMNAME) %>%
# { if (by_epu) dplyr::group_by(., EPU, add = T) else . } %>%
# { if (by_taxon) dplyr::group_by(., order, add = T) else . } %>%
# { if (by_season) dplyr::group_by(., SEASON, add = T) else . } %>%
dplyr::summarise(
# spawners_abund_lag0 = mean(spawners_abund,
# na.rm = T),
# spawners_abund_lag0_anom = mean(spawners_abund_lag0_anom,
# na.rm = T),
spawners_biom_lag0 = mean(spawners_biom, na.rm = T),
spawners_biom_lag0_anom = mean(spawners_biom_lag0_anom,
na.rm = T),
#recruits_abund = mean(recruits_abund, na.rm = T),
recruits_abund_lead1 = mean(recruits_abund_lead1,
na.rm = T),
recruits_abund_lead1_anom = mean(recruits_abund_lead1_anom,
na.rm = T),
# recruits_biom_lead1 = mean(recruits_biom_lead1,
# na.rm = T),
# recruits_biom_lead1_anom = mean(recruits_biom_lead1_anom,
# na.rm = T),
rs = mean(rs, na.rm = T),
# rs_abund = mean(rs_abund, na.rm = T),
# rs_biom = mean(rs_biom, na.rm = T),
rs_anom = mean(rs_anom, na.rm = T),
# rs_abund_anom = mean(rs_abund_anom, na.rm = T),
# rs_biom_anom = mean(rs_biom_anom, na.rm = T),
logr_abund_anom = mean(logr_abund_anom,
na.rm = T),
# logr_biom_anom = mean(logr_biom_anom,
# na.rm = T),
# logs_abund_anom = mean(logs_abund_anom,
# na.rm = T),
logs_biom_anom = mean(logs_biom_anom,
na.rm = T),
logrs_anom = mean(logrs_anom, na.rm = T),
# logrs_abund_anom = mean(logrs_abund_anom,
# na.rm = T),
# logrs_biom_anom = mean(logrs_biom_anom,
# na.rm = T),
# spawner_wl_anom = mean(spawner_wl_anom,
# na.rm = T),
# recruit_wl_lead1_anom = mean(recruit_wl_lead1_anom,
# na.rm = T),
# spawner_wl_res_anom = mean(spawner_wl_res_anom,
# na.rm = T),
# recruit_wl_res_lead1_anom = mean(recruit_wl_res_lead1_anom,
# na.rm = T)
)
saveRDS(dat_spec_rec, "AssessFishProdAnomaly.rds")Use Perretti code to make anomaly plots
Sub in interactive plots from here
Note that the code has been modified to only plot the black “total” line for years when all plotted stocks are assessed. At present this limits the total line for all stocks to the time series for GB haddock, 2010-2017. The mid atlantic stocks have more species with longer time series so the black line spans more of the plot.
Recruits/Spawner Biomass anomaly
# modify from 2-load,R
dat_spec_rec_SS <-
dat_spec_rec %>%
ungroup %>%
dplyr::select(YEAR, StockName, rs_anom) %>%
dplyr::rename(Time = YEAR,
Value = rs_anom) %>%
dplyr::mutate(Units = "anomaly (r/s)",
Var = StockName,
Source = "stocksmart") %>%
dplyr::select(-StockName)
#### Adjust plot properties -------------------------------
adjustAxes <-
ggplot2::theme(axis.title = element_text(size = 18),
axis.text = element_text(size = 15),
plot.title = element_text(size = 20))
### Plot function with interactive bars to see which species
#### Plot stacked bar with cpts for single var ------------
plot_stackbarcpts_single <- function(YEAR, var2bar,
x, xlab, ylab,
titl,
file_suffix,
leg_font_size = 10,
remove_leg = FALSE,
leg_ncol = 1,
wcpts = TRUE,
wdashed = TRUE,
height = 5.5,
width = 8,
filt = TRUE,
label = label,
y.text = y.text,
aggregate = FALSE,
intplot = FALSE) {
dat2bar <- data.frame(YEAR, var2bar,
x)
if (filt == TRUE){mab_species <- list("SUMMER FLOUNDER","SCUP","BLACK SEA BASS","BLUEFISH",
"NORTHERN SHORTFIN SQUID", "LONGFIN SQUID", "ATLANTIC MACKEREL",
"BUTTERFISH","ATLANTIC SURFCLAM", "OCEAN QUAHOG", "TILEFISH",
"BLUELINE TILEFISH","SPINY DOGFISH", "GOOSEFISH")
dat2plot <-
dat2bar %>%
tidyr::gather(variable, value, -YEAR, -var2bar) %>%
dplyr::mutate(var2bar = gsub(pattern = "_",
replacement = " ",
x = var2bar),
var2bar = gsub(pattern = "Atl.",
replacement = "ATLANTIC",
x = var2bar),
var2bar = gsub(pattern = "Atl",
replacement = "ATLANTIC",
x = var2bar),
var2bar = gsub(pattern = "NS and combined",
replacement = "",
x = var2bar),
var2bar = gsub(pattern = "YT",
replacement = "Yellowtail",
x = var2bar),
var2bar = gsub(pattern = " GoM",
replacement = " GOM",
x = var2bar),
var2bar = gsub(pattern = " by EPU",
replacement = "",
x = var2bar)) %>%
filter(var2bar %in% mab_species)
} else if (filt == FALSE){
dat2plot <-
dat2bar %>%
tidyr::gather(variable, value, -YEAR, -var2bar) %>%
dplyr::mutate(var2bar = gsub(pattern = "_",
replacement = " ",
# x = var2bar),
# var2bar = gsub(pattern = "Atl.",
# replacement = "ATLANTIC",
# x = var2bar),
# var2bar = gsub(pattern = "Atl",
# replacement = "ATLANTIC",
# x = var2bar),
# var2bar = gsub(pattern = "NS and combined",
# replacement = "",
# x = var2bar),
# var2bar = gsub(pattern = "YT",
# replacement = "Yellowtail",
# x = var2bar),
# var2bar = gsub(pattern = " GoM",
# replacement = " GOM",
# x = var2bar),
# var2bar = gsub(pattern = " by EPU",
# replacement = "",
x = var2bar))
}
if (aggregate){
agg <- dat2plot %>%
group_by(YEAR) %>%
dplyr::summarise(Total = sum(value, na.rm = T),
Count = n()) %>% # SG add a count of species
mutate(Totalold = ifelse(Total == 0, NA, Total),
Total = ifelse(Count < max(Count), NA, Total)) # SG keep only for max number species
}
p <-
ggplot(dat2plot,
aes(x = YEAR)) +
{if(intplot) geom_bar_interactive(data = dat2plot %>% filter(value > 0),
aes(y = value, fill = var2bar,
tooltip = var2bar, data_id = var2bar),
stat = "identity")} +
{if(intplot) geom_bar_interactive(data = dat2plot %>% filter(value < 0),
aes(y = value, fill = var2bar,
tooltip = var2bar, data_id = var2bar),
stat = "identity")} +
{if(!intplot) geom_bar(data = dat2plot %>% filter(value > 0),
aes(y = value, fill = var2bar),
stat = "identity")} +
{if(!intplot) geom_bar(data = dat2plot %>% filter(value < 0),
aes(y = value, fill = var2bar),
stat = "identity")} +
{if(aggregate) geom_line(data = agg,aes(x = YEAR, y = Total),
size = 1)} +
geom_hline(size = 0.3, aes(yintercept = 0)) +
xlab(xlab) +
ylab(ylab) +
ggtitle(titl) +
guides(fill = guide_legend(ncol = leg_ncol)) +
theme_ts()+
theme(axis.title = element_text(size = 16),
axis.text = element_text(size = 15),
plot.title = element_text(size = 20),
legend.text = element_text(size = leg_font_size),
legend.title = element_blank()) +
annotate("text", label = label, x = 1980, y = y.text,size = 8, colour = "black")
if(remove_leg) p <- p + theme(legend.position = "none")
return(p)
}
# bar_dat <- ecodata::productivity_anomaly %>%
# filter(EPU == "MAB")
# mafmc <-plot_stackbarcpts_single(YEAR = bar_dat$Time,
# var2bar = bar_dat$Var,
# x = bar_dat$Value,
# titl = "",
# xlab = "",
# ylab = "Small fish per large fish biomass (anomaly)",
# height = 5.5,
# width = 9,
# filt = TRUE,
# label = "A",
# y.text = 4.5)
#
bar_dat <- dat_spec_rec_SS
allstocks <- plot_stackbarcpts_single(YEAR = bar_dat$Time,
var2bar = bar_dat$Var,
x = bar_dat$Value,
titl = "All stocks",
xlab = "",
ylab = "Recruits/Mature biomass (anomaly)",
height = 5.5,
width = 9,
filt = FALSE,
label = "",
y.text = 10,
aggregate = TRUE,
intplot = TRUE,
remove_leg = TRUE)
#all
ggiraph(code = print(allstocks), pointsize = 14, width_svg = 8, height_svg = 6)bar_dat <- dat_spec_rec_SS %>%
filter(Var %in% c("Atlantic mackerel - Gulf of Maine / Cape Hatteras",
"Tilefish - Mid-Atlantic Coast" ,
"Yellowtail flounder - Southern New England / Mid-Atlantic",
"Winter flounder - Southern New England / Mid-Atlantic" ,
"Summer flounder - Mid-Atlantic Coast" ,
"Scup - Atlantic Coast",
"Bluefish - Atlantic Coast",
"Black sea bass - Mid-Atlantic Coast",
"Butterfish - Gulf of Maine / Cape Hatteras",
"Spiny dogfish - Atlantic Coast"
))
midstocks <- plot_stackbarcpts_single(YEAR = bar_dat$Time,
var2bar = bar_dat$Var,
x = bar_dat$Value,
titl = "Mid-Atlantic stocks",
xlab = "",
ylab = "Recruits/Mature biomass (anomaly)",
height = 5.5,
width = 9,
filt = FALSE,
label = "",
y.text = 10,
aggregate = TRUE,
intplot = TRUE,
remove_leg = TRUE)
#all
ggiraph(code = print(midstocks), pointsize = 14, width_svg = 8, height_svg = 6)For Mid-Atlantic stocks plot I am including:
- “Atlantic mackerel - Gulf of Maine / Cape Hatteras”
- “Tilefish - Mid-Atlantic Coast”
- “Yellowtail flounder - Southern New England / Mid-Atlantic”
- “Winter flounder - Southern New England / Mid-Atlantic”
- “Summer flounder - Mid-Atlantic Coast”
- “Scup - Atlantic Coast”
- “Bluefish - Atlantic Coast”
- “Black sea bass - Mid-Atlantic Coast”
- “Butterfish - Gulf of Maine / Cape Hatteras”
- “Spiny dogfish - Atlantic Coast”
Below the static plot code, still needs work
# static plot, needs legend tweakes
# plot functions
#source('~/Documents/0_Data/ESR/SOE2022/trawlr/run_for_soe/1-func.R')
# modify from 2-load,R
# dat_spec_rec_SS <-
# dat_spec_rec %>%
# ungroup %>%
# dplyr::select(YEAR, StockName, rs_anom) %>%
# dplyr::rename(Time = YEAR,
# Value = rs_anom) %>%
# dplyr::mutate(Units = "anomaly (r/s)",
# Var = StockName,
# Source = "stocksmart") %>%
# dplyr::select(-StockName)
#
# # modify from 4-do.R
# ## Recruit per spawner (all stocks in one figure) --------
# plot_stackbarcpts_single(YEAR = dat_spec_rec_SS$Time,
# var2bar = dat_spec_rec_SS$Var,
# x = dat_spec_rec_SS$Value,
# titl = "",
# xlab = "",
# ylab = "Recruits/Mature biomass (anomaly)",
# height = 5.5,
# width = 9)Recruits anomaly
dat_spec_reconly_SS <-
dat_spec_rec %>%
ungroup %>%
dplyr::select(YEAR, StockName, recruits_abund_lead1_anom) %>%
dplyr::rename(Time = YEAR,
Value = recruits_abund_lead1_anom) %>%
dplyr::mutate(Units = "anomaly",
Var = StockName,
Source = "stocksmart") %>%
dplyr::select(-StockName)
bar_dat <- dat_spec_reconly_SS
allstocks <- plot_stackbarcpts_single(YEAR = bar_dat$Time,
var2bar = bar_dat$Var,
x = bar_dat$Value,
titl = "All stocks",
xlab = "",
ylab = "Recruits (anomaly)",
height = 5.5,
width = 9,
filt = FALSE,
label = "",
y.text = 10,
aggregate = TRUE,
intplot = TRUE,
remove_leg = TRUE)
#all
ggiraph(code = print(allstocks), pointsize = 14, width_svg = 8, height_svg = 6)bar_dat <- dat_spec_reconly_SS %>%
filter(Var %in% c("Atlantic mackerel - Gulf of Maine / Cape Hatteras",
"Tilefish - Mid-Atlantic Coast" ,
"Yellowtail flounder - Southern New England / Mid-Atlantic",
"Winter flounder - Southern New England / Mid-Atlantic" ,
"Summer flounder - Mid-Atlantic Coast" ,
"Scup - Atlantic Coast",
"Bluefish - Atlantic Coast",
"Black sea bass - Mid-Atlantic Coast",
"Butterfish - Gulf of Maine / Cape Hatteras"
))
midstocks <- plot_stackbarcpts_single(YEAR = bar_dat$Time,
var2bar = bar_dat$Var,
x = bar_dat$Value,
titl = "Mid-Atlantic stocks",
xlab = "",
ylab = "Recruits (anomaly)",
height = 5.5,
width = 9,
filt = FALSE,
label = "",
y.text = 10,
aggregate = TRUE,
intplot = TRUE,
remove_leg = TRUE)
#all
ggiraph(code = print(midstocks), pointsize = 14, width_svg = 8, height_svg = 6)Spawner Biomass anomaly
dat_spec_ssbonly_SS <-
dat_spec_rec %>%
ungroup %>%
dplyr::select(YEAR, StockName, spawners_biom_lag0_anom) %>%
dplyr::rename(Time = YEAR,
Value = spawners_biom_lag0_anom) %>%
dplyr::mutate(Units = "anomaly",
Var = StockName,
Source = "stocksmart") %>%
dplyr::select(-StockName)
bar_dat <- dat_spec_ssbonly_SS
allstocks <- plot_stackbarcpts_single(YEAR = bar_dat$Time,
var2bar = bar_dat$Var,
x = bar_dat$Value,
titl = "All stocks",
xlab = "",
ylab = "Mature biomass (anomaly)",
height = 5.5,
width = 9,
filt = FALSE,
label = "",
y.text = 10,
aggregate = TRUE,
intplot = TRUE,
remove_leg = TRUE)
#all
ggiraph(code = print(allstocks), pointsize = 14, width_svg = 8, height_svg = 6)bar_dat <- dat_spec_ssbonly_SS %>%
filter(Var %in% c("Atlantic mackerel - Gulf of Maine / Cape Hatteras",
"Tilefish - Mid-Atlantic Coast" ,
"Yellowtail flounder - Southern New England / Mid-Atlantic",
"Winter flounder - Southern New England / Mid-Atlantic" ,
"Summer flounder - Mid-Atlantic Coast" ,
"Scup - Atlantic Coast",
"Bluefish - Atlantic Coast",
"Black sea bass - Mid-Atlantic Coast",
"Butterfish - Gulf of Maine / Cape Hatteras"
))
midstocks <- plot_stackbarcpts_single(YEAR = bar_dat$Time,
var2bar = bar_dat$Var,
x = bar_dat$Value,
titl = "Mid-Atlantic stocks",
xlab = "",
ylab = "Mature biomass (anomaly)",
height = 5.5,
width = 9,
filt = FALSE,
label = "",
y.text = 10,
aggregate = TRUE,
intplot = TRUE,
remove_leg = TRUE)
#all
ggiraph(code = print(midstocks), pointsize = 14, width_svg = 8, height_svg = 6)References
Perretti, Ct, Mj Fogarty, Kd Friedland, Ja Hare, Sm Lucey, Rs McBride,
Tj Miller, et al. 2017. “Regime Shifts in Fish Recruitment on the
Northeast US Continental
Shelf.” Marine Ecology Progress Series 574
(July): 1–11. https://doi.org/10.3354/meps12183.