Get SOE stock status indicator from assessment data

The new way

Andy made the assessmentdata package based on Stock SMART. Surely this will be better?

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

In stockAssessmentData we have time series. Columns are Species, Region, Year, Value, Metric, Description, Units, AssessmentYear and the reported metrics are Catch, Fmort, Recruitment, Abundance, Index.

datatable(head(stockAssessmentData), rownames = FALSE)

In stockAssessmentSummary we have assessment metadata. Columns are Stock Name, Jurisdiction, FMP, Science Center, Regional Ecosystem, FSSI Stock?, ITIS Taxon Serial Number, Scientific Name, Common Name, Stock Area, Assessment Year, Assessment Month, Last Data Year, Update Type, 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, Assessment Type, 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.

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

Trying to rebuild my spreadsheet from stockAssessmentSummary and use the ecodata code to make the dataset:

test2019assess <- 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`)

write.csv(test2019assess, here("stockSMART2019assess.csv"))

# from get_stocks.R, ecodata 2020
  #assess <- read.csv(file.path(data.dir, "2019assess.csv"))
  assess <- test2019assess
  #decode <- read.csv(file.path(data.dir, "2019decoder.csv"))
  decode <- read.csv(here("2019decoder.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 I get the same plot for assessments up to Dec 2019 I got the old way.

Comparisons

StockSMART 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.0 #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"),
                    x = rep(0.9*x.max,5), y = seq(0.93*y.max,1.4,-.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 2020 ecodata source, Mid-Atlantic

stock_status <- 
  ecodata::stock_status %>%
  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.0 #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"),
                    x = rep(0.9*x.max,5), y = seq(0.93*y.max,1.4,-.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()

StockSMART 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 = 8))
b.df <- data.frame(text = b_missing$Code,
                    x = rep(x.max*0.7,length(b_missing$Code)),
                    y = c(1.45,2.15))
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 = 8))

#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 2020 ecodata source, New England

stock_status <- ecodata::stock_status %>%
  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))
b.df <- data.frame(text = b_missing$Code,
                    x = rep(x.max*0.7,length(b_missing$Code)),
                    y = c(1.45,2.15))
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 = 8))

#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()

Issues

Mid-Atlantic: Mackerel B.Bmsy is different, lower in stockSMART. New England: GOM Cod called F and B missing in stockSMART. Hakes are different. Wolffish missing from stockSMART?

Comparing the 2019assess.csv with the stockSMART2019assess.csv, discrepancies are:

Mackerel is correct in stockSMART. The 2019assess.csv spreadsheet has Blimit entered where Bmsy should be.

Next…

Then apply for assessments up to Dec 2020. Well, max assessment year is 2019, so I guess I still need to fill this in by hand!

Options:

Retain only the fields used by get_stocks.R so maybe I don’t break ecodata.
Recreate the entire spreadsheet for all assessments (means they must all be in stocksmart).

Decisions

Update with stockSMART for 2019 and prior records, add only the 2020 assessments by hand. (Is stockSMART the definitive source of status information?)
Recreate assess spreadsheet fields from stockSMART including all fields to act as metadata

new2019assess <- stockAssessmentSummary %>%
  filter(`Science Center` == "NEFSC") %>%
  rename(Entity.Name = `Stock Name`) %>%
  rename_all(list(~make.names(.)))

Add 2020 assessments using a script here rather than typing into a spreadsheet

#subset stockSMARTsummary for updates species
  #add a row for each species with an update, copying most recent data
  #update Assessment.Year = 2020
  #update all other fields--nope just got them already compiled!
  
# updatespecies <- new2019assess %>%
#   filter(str_detect(Entity.Name, c("redfish|halibut|herring|scallop|surfclam|wolffish|Butterfish|Longfin|pout|quahog|Red hake|Silver hake|Windowpane|Winter flounder"))) %>%
#   group_by(Entity.Name) %>%
#   do(add_row(., Assessment.Year = 2020)) %>%
#   ungroup() %>%
#   fill(!c(Assessment.Year, Assessment.Month, Last.Data.Year,
#           Update.Type, Review.Result, Assessment.Model, Model.Version,
#           Citation:MSY.Unit))

prelim2020 <- read.csv(here("NE Stock Assessment Results.csv")) %>%
  filter(Assessment.Year == 2020) %>%
  rename(Entity.Name = Stock,
         FSSI.Stock. = FSSI,
         Estimated.F = Best.F,
         Estimated.B = Best.B,
         Review.Result = Review.Type) %>%
  select(-c(Year, Status.Stock., Record.Status, TimeSeries.Data.,
            Survey.Links., Adequate, Minimum.F, Maximum.F,
            Minimum.B, Maximum.B, 
            Stock.Level.Relative.to.Bmsy:Decision.memo.related.to.inadequate.rebuilding.progress))


update2020assess <- bind_rows(new2019assess, prelim2020)

write.csv(update2020assess, here("2020assess.csv"))

Retain only Entity.Name, Council, and Code fields (used by get_stocks) in decoder

newdecoder <- read.csv(here("2019decoder.csv")) %>%
  select(Entity.Name, Code, Council)

write.csv(newdecoder, here("2020decoder.csv"))

2020 stock assessment source

Used the form here: https://apps-nefsc.fisheries.noaa.gov/saw/sasi/sasi_report_options.php

Files should be up-to-date as of the day of my query, December 9, 2020, according to the website.

Checked the boxes: Year–2020

Checked each available 2020 species and stock area in turn, downloaded .zip of “all files”.

Species with 2020 updates: Acadian redfish, Atlantic halibut, Atlantic herring, Atlantic Sea Scallop, Atlantic surfclam, Atlantic wolffish, Butterfish, Longfin squid, Ocean Pout, Ocean quahog, Red Hake (2 stocks), Silver hake (2 stocks), Windowpane flounder (2 stocks), Winter flounder (3 stocks).

Update: these results were compiled as preliminary information by Jeff Vieser, who is awaiting final NEFSC publication to post to stockSMART. He sent me the spreadsheet NE Stock Assessment Results.xlsx 10 December 2020. We will mark these preliminary, pending final number confirmation.

I will spot check the spreadsheet to ensure it matches numbers in the downloaded reports.

New Plots with preliminary 2020 assessment results

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

stock_status_2020 <-
    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))

StockSMART source, Mid-Atlantic

stock_status <- 
  stock_status_2020 %>%
  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.0 #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"),
                    x = rep(0.9*x.max,5), y = seq(0.93*y.max,1.4,-.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()

StockSMART source, New England

stock_status <- 
  stock_status_2020 %>%
  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.0, 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.35, 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.2, length.out = length(f_missing$Code)))

#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()

Get SOE stock status indicator from assessment data

Sarah Gaichas

10 December, 2020

The old way

The new way

Comparisons

StockSMART source, Mid-Atlantic

SOE 2020 ecodata source, Mid-Atlantic

StockSMART source, New England

SOE 2020 ecodata source, New England

Issues

Next…

Decisions

2020 stock assessment source

New Plots with preliminary 2020 assessment results

StockSMART source, Mid-Atlantic

StockSMART source, New England