Major Takeaway

The stockSMART updated plots and my by-hand dataset for 2023 SOE are the same except I included the status results from the December 2022 research track assessments for spiny dogfish and bluefish, both of which changed status. We need to decide whether we want these in the SOE or not.

Pull from stocksmart for 2024 SOE reports

Andy renamed the assessmentdata package stocksmart based on Stock SMART.

Two data frames are in the package, stockAssessmentData and stockAssessmentSummary.

In stockAssessmentData we have time series. Columns are StockName, Stockid, Assessmentid, Year, Value, Metric, Description, Units, AssessmentYear, Jurisdiction, FMP, CommonName, ScientificName, ITIS, AssessmentType, StockArea, RegionalEcosystem and the reported metrics are Catch, Fmort, Recruitment, Abundance, Index.

datatable(head(stockAssessmentData), rownames = FALSE)

In stockAssessmentSummary we have assessment metadata. Columns are Stock ID, Stock Name, Jurisdiction, FMP, Science Center, Regional Ecosystem, FSSI Stock?, ITIS Taxon Serial Number, Scientific Name, Common Name, Stock Area, Assessment ID, Assessment Year, Assessment Month, Last Data Year, Review Result, Assessment Model, Model Version, Lead Lab, Citation, Final Assessment Report 1, Final Assessment Report 2, Point of Contact, Life History Data, Abundance Data, Catch Data, Assessment Level, Assessment Frequency, Model Category, Catch Input Data, Abundance Input Data, Biological Input Data, Ecosystem Linkage, Composition Input Data, F Year, Estimated F, F Unit, F Basis, Flimit, Flimit Basis, Fmsy, Fmsy Basis, F/Flimit, F/Fmsy, Ftarget, Ftarget Basis, F/Ftarget, B Year, Estimated B, B Unit, B Basis, Blimit, Blimit Basis, Bmsy, Bmsy Basis, B/Blimit, B/Bmsy, MSY, MSY Unit, Assessment Type.

datatable(head(stockAssessmentSummary), rownames = FALSE, options = list(scrollX = TRUE))

Build ecodata input spreadsheet from stockAssessmentSummary and use the ecodata code to make the dataset for plotting:

assess2023 <- stockAssessmentSummary %>%
  filter(`Science Center` == "NEFSC") %>%
  select(c(`Stock Name`, Jurisdiction, FMP, `Science Center`, 
           `Stock Area`, `Assessment Year`, `Last Data Year`,
           `F Year`, `Estimated F`, Flimit, Fmsy, `F/Flimit`, 
           `F/Fmsy`, Ftarget, `F/Ftarget`, `B Year`, `Estimated B`,
           `B Unit`, Blimit, Bmsy, `B/Blimit`, `B/Bmsy`)) %>%
  arrange(Jurisdiction, `Stock Name`, FMP, `Assessment Year`) %>%
  rename(Entity.Name = `Stock Name`,
         Assessment.Year = `Assessment Year`,
         F.Fmsy = `F/Fmsy`,
         B.Bmsy = `B/Bmsy`,
         Estimated.F = `Estimated F`,
         Estimated.B = `Estimated B`)

write.csv(assess2023, here("assess.csv"))

# from get_stocks.R, ecodata 2020
  #assess <- read.csv(file.path(data.dir, "2019assess.csv"))
  assess <- assess2023
  #decode <- read.csv(file.path(data.dir, "2019decoder.csv"))
  decode <- read.csv(here("2020decoder.csv"))
  
write.csv(decode, here("decoder.csv"))

stock_status_stockSMART <-
    assess %>%
    dplyr::group_by(Entity.Name) %>%
    dplyr::filter(Assessment.Year == max(Assessment.Year)) %>%
  #Find last year assessment occurred for each stock
    dplyr::ungroup() %>%
    dplyr::left_join(.,decode, by = "Entity.Name") %>% #Join in list of managed species
    dplyr::select(Entity.Name, Assessment.Year, F.Fmsy, B.Bmsy, Council, Code) %>%
  #select column variables to keep
    dplyr::mutate(id = 1:length(Entity.Name)) %>%
    tidyr::gather(.,Var, Value,-id,-Entity.Name,-Assessment.Year,-Council,-Code) %>%
  #wide to long
    dplyr::select(-id) %>%
    dplyr::rename(`Last assessment` = Assessment.Year,
                  Stock = Entity.Name) %>% #rename variables for clarity
    dplyr::mutate(Units = "unitless") #%>%
    #dplyr::mutate(Value = replace(Value, which(Code == "N Windowpane" & Var == "F.Fmsy"), NA))

Then test to see if we see the updates relative to SOE 2023 (updates through 2022) ecodata. I’m leaving out all the plot annotations for unknown status.

Comparisons

StockSMART Oct 2023 source, Mid-Atlantic

stock_status <- 
  stock_status_stockSMART %>%
  mutate(Code = recode(Code, "Dogfish" = "Sp. Dogfish" )) %>% 
  spread(.,Var,Value) %>% 
  filter(Council %in% c("MAFMC","Both")) %>% 
  group_by(Stock) %>% 
  mutate(score = case_when(
    (B.Bmsy <0.5) ~"a",
    (F.Fmsy >1) ~ "a", 
    (F.Fmsy < 1 & B.Bmsy > 0.5 & B.Bmsy < 1) ~ "b",
    (F.Fmsy < 1 & B.Bmsy > 1) ~ "c"))
#Plot constants
y.max <- 2.1 #2.0 mackerel cut off F/Fmsy is 2.08
x.max <- 2.6
#A dataframe that defines custom legend for stocks with unknown status
# unknown <- data.frame(text = c("Unknown Status", "Longfin Squid",
#                               "Shortfin Squid", "N. Goosefish", "S. Goosefish"),
#                     x = rep(0.9*x.max,5), y = seq(0.93*y.max,1.5,-.1))

# Custom Color
custom_color<- c("#56B4E9", "#009E73", "#0072B2")
#Plotting code
ggplot(data = stock_status) +
  geom_vline(xintercept = 1, linetype = "dotted")+
  geom_vline(xintercept = 0.5, linetype = "dashed")+
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_point(aes(x = B.Bmsy,
                 y = F.Fmsy,
                 shape = Council,
                 color = score)) +
  geom_text_repel(aes(x = B.Bmsy, #geom_text_repel auto-jitters text around points
                      y = F.Fmsy,
                      label = Code, 
                      color = score), 
                  show.legend = FALSE, nudge_y = -0.01, nudge_x = 0.05) +
  scale_color_brewer(palette = "Dark2",
                     breaks = stock_status$score) +
  ylim(0,y.max) +
  xlim(0,x.max) +
  # geom_text(data = unknown, aes(x = x, y = y, label = text), #Custom legend for unknown stock status
  #           size = c(4.75,rep(4,4))) +
  # annotate("rect", xmin = 0.8*x.max,
  #          xmax = x.max,
  #          ymin = 0.65*y.max,
  #          ymax = 0.90*y.max,
  #          alpha = 0.1) +
  xlab(expression(~B/B[msy])) +
  ylab(expression(~F/F[msy])) +
  guides(color = FALSE) +
  theme_ts()

SOE 2023 ecodata source, Mid-Atlantic

#Get data, spread for plotting, and filter
stock_status <- ecodata::stock_status %>%
  dplyr::mutate(Code = recode(Code, "Dogfish" = "Sp. Dogfish" ), 
                Code = recode(Code, "Mackerel" = "At. Mackerel")) %>% 
  tidyr::spread(.,Var,Value) %>% 
  dplyr::filter(Council %in% c("MAFMC","Both")) %>% 
  dplyr::group_by(Stock) %>% 
  dplyr::mutate(score = case_when(
    (B.Bmsy <0.5) ~"a",
    (F.Fmsy >1) ~ "a", 
    (F.Fmsy < 1 & B.Bmsy > 0.5 & B.Bmsy < 1) ~ "b",
    (F.Fmsy < 1 & B.Bmsy > 1) ~ "c"))
#Plot constants
y.max <- 2.1 #1.75 mackerel cut off F/Fmsy is 1.8
x.max <- 2.6
#A dataframe that defines custom legend for stocks with unknown status
unknown <- data.frame(text = c("Unknown Status", "Longfin Squid",
                              "Shortfin Squid", "N. Goosefish", "S. Goosefish", "Blueline Tilefish", "Chub Mackerel"),
                    x = rep(0.9*x.max,7), y = seq(0.88*y.max,1.2,-0.1))

# Custom Color
custom_color<- c("#56B4E9", "#009E73", "#0072B2")
#Plotting code
ggplot2::ggplot(data = stock_status) +
  ggplot2::geom_vline(xintercept = 1, linetype = "dotted")+
  ggplot2::geom_vline(xintercept = 0.5, linetype = "dashed")+
  ggplot2::geom_hline(yintercept = 1, linetype = "dashed") +
  ggplot2::geom_point(aes(x = B.Bmsy,
                 y = F.Fmsy,
                 shape = Council,
                 color = score)) +
  ggrepel::geom_text_repel(aes(x = B.Bmsy, #geom_text_repel auto-jitters text around points
                      y = F.Fmsy,
                      label = Code, 
                      color = score), 
                  show.legend = FALSE, nudge_y = -0.01, nudge_x = 0.05) +
  ggplot2::scale_color_brewer(palette = "Dark2",
                     breaks = stock_status$score) +
  ggplot2::ylim(0,y.max) +
  ggplot2::xlim(0,x.max) +
  ggplot2::geom_text(data = unknown, aes(x = x-0.5, y = y+0.2, label = text), #Custom legend for unknown stock status
            size = c(4.75,rep(4,6))) +
  #ggplot2::annotate("rect", xmin = 0.8*x.max,
  #         xmax = x.max
  #         ymin = 0.65*y.max,
  #         ymax = 0.90*y.max,
  #         alpha = 0.1) +
  ggplot2::xlab(expression(~B/B[msy])) +
  ggplot2::ylab(expression(~F/F[msy])) +
  ggplot2::guides(color = FALSE) +
  ecodata::theme_ts()

StockSMART Oct 2023 source, New England

stock_status <- stock_status_stockSMART %>%
  mutate(Code = recode(Code, "Dogfish" = "Sp. Dogfish" )) %>% 
  spread(.,Var,Value) %>% 
  filter(Council %in% c("NEFMC","Both")) %>% 
  group_by(Stock) %>% 
  mutate(score = case_when(
    (B.Bmsy <0.5) ~"a",
    (F.Fmsy >1) ~ "a", 
    (F.Fmsy < 1 & B.Bmsy > 0.5 & B.Bmsy < 1) ~ "b",
    (F.Fmsy < 1 & B.Bmsy > 1) ~ "c"))
#Plot constants
y.max <- 1.5
x.max <- 10
all_missing <- stock_status %>%
  filter(is.na(B.Bmsy),is.na(F.Fmsy)) %>% 
  dplyr::select(Code, Council)
b_missing <- stock_status %>%
  filter(is.na(B.Bmsy), !is.na(F.Fmsy)) %>% 
  dplyr::select(Code, Council)
f_missing <- stock_status %>%
  filter(is.na(F.Fmsy), !is.na(B.Bmsy)) %>% 
  dplyr::select(Code, Council)
#A dataframe that defines custom legend for stocks with unknown status
# all.df <- data.frame(text = all_missing$Code,
#                     x = rep(x.max*0.9,length(all_missing$Code)),
#                     #y = seq(1.45,1.05, length.out = 7))
#                     y = seq(1.45,1.05, length.out = length(all_missing$Code)))
# b.df <- data.frame(text = b_missing$Code,
#                     x = rep(x.max*0.7,length(b_missing$Code)),
#                     y = c(1.45,2.15, length.out = length(b_missing$Code)))
# f.df <- data.frame(text = f_missing$Code,
#                     x = rep(x.max*0.5,length(f_missing$Code)),
#                     y = seq(1.45,1.0, length.out = length(f_missing$Code)))

# Custom Color
custom_color<- c("#56B4E9", "#009E73", "#0072B2")

#Plotting code
ggplot(data = stock_status) +
  geom_vline(xintercept = 1, linetype = "dotted", color = "grey60")+
  geom_vline(xintercept = 0.5, linetype = "dashed", color = "grey60")+
  geom_hline(yintercept = 1, linetype = "dashed", color = "grey60") +
  geom_point(aes(x = B.Bmsy,
                 y = F.Fmsy,
                 color = stock_status$score)) +
  geom_text_repel(aes(x = B.Bmsy, #geom_text_repel auto-jitters text around points
                      y = F.Fmsy,
                      label = Code,
                      color = stock_status$score), show.legend = FALSE,nudge_y = -0.01, nudge_x = 0.05) +
  ylim(0,y.max) +
  xlim(0,x.max*1.1) +
  # geom_text(data = all.df, aes(x = x, y = y, label = text),show.legend = FALSE, size = 3)+
  # geom_text(data = b.df, aes(x = x, y = y, label = text),show.legend = FALSE, size = 3)+
  # geom_text(data = f.df, aes(x = x, y = y, label = text),show.legend = FALSE, size = 3)+
  # scale_color_brewer(palette = "Dark2", #Change legend labels for clarity
  #                    breaks = stock_status$score) +
  # annotate("rect", xmin = 0.924*x.max,
  #          xmax = 1.08*x.max,
  #          ymin = 0.645*y.max,
  #          ymax = 0.98*y.max,
  #          alpha = 0.01) +
  # annotate("text", x = 9, y = 1.5, label = "F and B missing", fontface =2, size = 3)+
  # annotate("rect",  
  #            xmin = 0.70*x.max,
  #            xmax = 0.85*x.max,
  #            ymin = 0.90*y.max,
  #            ymax = 1.8,
  #            alpha = 0.01) +
  # annotate("text", x = 7, y = 1.5, label = "B missing", fontface =2, size = 3)+
  # annotate("rect", xmin = 0.509*x.max,
  #          xmax = 0.681*x.max,
  #          ymin = 0.65*y.max,
  #          ymax = 0.98*y.max,
  #          alpha = 0.01) +
  # annotate("text", x = 5, y = 1.5, label = "F missing", fontface =2, size = 3)+
  xlab(expression(~B/B[msy])) +
  ylab(expression(~F/F[msy])) +
  guides(color = FALSE) +
  theme_ts()

SOE 2023 ecodata source, New England

#Get data, spread for plotting, and filter
stock_status <- ecodata::stock_status %>%
  dplyr::mutate(Code = dplyr::recode(Code, "Dogfish" = "Sp. Dogfish" )) %>% 
  tidyr::spread(.,Var,Value) %>% 
  dplyr::filter(Council %in% c("NEFMC","Both")) %>% 
  dplyr::group_by(Stock) %>% 
  dplyr::mutate(score = case_when(
    (B.Bmsy <0.5) ~"a",
    (F.Fmsy >1) ~ "a", 
    (F.Fmsy < 1 & B.Bmsy > 0.5 & B.Bmsy < 1) ~ "b",
    (F.Fmsy < 1 & B.Bmsy > 1) ~ "c"))
#Plot constants
y.max <- 1.5
x.max <- 10
all_missing <- stock_status %>%
  dplyr::filter(is.na(B.Bmsy),is.na(F.Fmsy)) %>% 
  dplyr::select(Code, Council)
b_missing <- stock_status %>%
  dplyr::filter(is.na(B.Bmsy), !is.na(F.Fmsy)) %>% 
  dplyr::select(Code, Council)
f_missing <- stock_status %>%
  dplyr::filter(is.na(F.Fmsy), !is.na(B.Bmsy)) %>% 
  dplyr::select(Code, Council)
#A dataframe that defines custom legend for stocks with unknown status
all.df <- data.frame(text = all_missing$Code,
                    x = rep(x.max*0.9,length(all_missing$Code)),
                    y = seq(1.45,0.7, length.out = length(all_missing$Code)))
b.df <- data.frame(text = b_missing$Code,
                    x = rep(x.max*0.7,length(b_missing$Code)),
                    y = seq(1.45,1.30, length.out = length(b_missing$Code)))
f.df <- data.frame(text = f_missing$Code,
                    x = rep(x.max*0.5,length(f_missing$Code)),
                    y = seq(1.45,1, length.out = length(f_missing$Code)))

#Plotting code
ggplot2::ggplot(data = stock_status) +
  ggplot2::geom_vline(xintercept = 1, linetype = "dotted", color = "grey60")+
  ggplot2::geom_vline(xintercept = 0.5, linetype = "dashed", color = "grey60")+
  ggplot2::geom_hline(yintercept = 1, linetype = "dashed", color = "grey60") +
  ggplot2::geom_point(aes(x = B.Bmsy,
                 y = F.Fmsy,
                 color = stock_status$score)) +
  ggrepel::geom_text_repel(aes(x = B.Bmsy, #geom_text_repel auto-jitters text around points
                      y = F.Fmsy,
                      label = Code,
                      color = stock_status$score), show.legend = FALSE,nudge_y = -0.01, nudge_x = 0.05) +
  ggplot2::ylim(0,y.max) +
  ggplot2::xlim(0,x.max*1.1) +
  ggplot2::geom_text(data = all.df, aes(x = x, y = y, label = text),show.legend = FALSE, size = 3)+
  ggplot2::geom_text(data = b.df, aes(x = x, y = y, label = text),show.legend = FALSE, size = 3)+
  ggplot2::geom_text(data = f.df, aes(x = x, y = y, label = text),show.legend = FALSE, size = 3)+
  ggplot2::scale_color_brewer(palette = "Dark2", #Change legend labels for clarity
                     breaks = stock_status$score) +
  ggplot2::annotate("rect", xmin = 0.924*x.max,
           xmax = 1.08*x.max,
           ymin = 0.645*y.max,
           ymax = 0.98*y.max,
           alpha = 0.01) +
  ggplot2::annotate("text", x = 9, y = 1.5, label = "F and B missing", fontface =2, size = 3)+
  ggplot2::annotate("rect",  
             xmin = 0.70*x.max,
             xmax = 0.85*x.max,
             ymin = 0.30*y.max,
             ymax = 0.98*y.max,
             alpha = 0.01) +
  ggplot2::annotate("text", x = 7, y = 1.5, label = "B missing", fontface =2, size = 3)+
  ggplot2::annotate("rect", xmin = 0.509*x.max,
           xmax = 0.681*x.max,
           ymin = 0.65*y.max,
           ymax = 0.98*y.max,
           alpha = 0.01) +
  ggplot2::annotate("text", x = 5, y = 1.5, label = "F missing", fontface =2, size = 3)+
  ggplot2::xlab(expression(~B/B[msy])) +
  ggplot2::ylab(expression(~F/F[msy])) +
  ggplot2::guides(color = FALSE) +
  ecodata::theme_ts()

Issues

None noted, stockSMART and stocksmart have been updated as of December 13 2023.