cjfas-2023-0093suppld.Rmd

---
title: "Assessing small pelagic fish trends in space and time using piscivore stomach contents: Supplement d, Sensitivities"
author: "Sarah Gaichas, James Gartland, Brian Smith, Anthony Wood, Elizabeth Ng, Michael Celestino,
  Katie Drew, Abigail Tyrell, and James Thorson"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  bookdown::word_document2:
    toc: true
  bookdown::html_document2: 
    toc: true
    toc_float: true
  bookdown::pdf_document2:
    includes: 
       in_header: latex/header1.tex
    keep_tex: true
link-citations: yes
csl: "canadian-journal-of-fisheries-and-aquatic-sciences.csl"
bibliography: FishDiet_EcoIndicators.bib
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, # no code blocks in word doc
                      message = FALSE,
                      warning = FALSE)

library(tidyverse)
theme_set(theme_bw())
library(here)

library(officedown)
library(officer)
library(flextable)
op_section <- prop_section(type = "continuous")
close_section <- prop_section(
    page_size = page_size(orient = "landscape"), 
    type = "continuous")

```

# Forage index sensitivity to included prey species {-}

## Different prey cut-off {-}

Two methods were applied to derive the bluefish prey list. Initially, a cut-off of 10 observations of a prey item across all sampled bluefish stomachs was applied separately to the NEFSC food habits data and the NEAMAP food habits data. This approach resulted in the same common bluefish prey items for each dataset, but slightly different uncommon prey items. For example, Atlantic herring, Atlantic mackerel, and silver hake were included as prey items for NEFSC data, but not for NEAMAP data using this approach, due to the timing of the surveys combined with the distribution of large and small bluefish and these prey items inshore versus offshore. 

The second method (the one used for this paper) first combined the prey observations across NEFSC and NEAMAP food habits datasets and then applied a cutoff of 25 observations of a prey item to the combined dataset. This ensured that a consistent set of prey were used across predators for both datasets. The advantage of this approach is that any Atlantic herring, Atlantic mackerel, and silver hake seen in non-bluefish predator stomachs in the NEAMAP dataset would still be included in the prey index to reflect these prey nearshore. Further, categories reflecting different survey naming conventions for the same prey were combined using this approach (e.g., *Illex illecebrosus* and *Illex* sp.)

While the second method used in this paper is appealing in its consistent treatment of prey items across the two datasets, the resulting index showed extremely low sensitivity to changing the prey cutoff (Figs. Sd\@ref(fig:preycut1)-Sd\@ref(fig:preycut2)). This suggests that the index is primarily driven by the most commonly identified prey items across the datasets.  

```{r preycut1, fig.cap="Trend comparison between fall forage indices using different bluefish prey cut-offs: 10 or more observations by survey and 25 observations combined across surveys (current method).", fig.height=8}
# compare 2022 prey list to paper update preylist for indices

stratlook <- data.frame(Stratum = c("Stratum_1",
                                    "Stratum_2",
                                    "Stratum_3",
                                    "Stratum_4",
                                    "Stratum_5",
                                    "Stratum_6",
                                    "Stratum_7",
                                    "Stratum_8",
                                    "Stratum_9",
                                    "Stratum_10",
                                    "Stratum_11",
                                    "Stratum_12",
                                    "Stratum_13",
                                    "Stratum_14",
                                    "Stratum_15"),
                        Region  = c("AllEPU", 
                                    "MABGB", 
                                    "MABGBstate", 
                                    "MABGBfed", 
                                    "MAB",
                                    "GB",
                                    "GOM",
                                    "bfall",
                                    "bfin",
                                    "bfoff",
                                    "MABGBalbinshore",
                                    "MABGBothoffshore",
                                    "albbfin",
                                    "albbfall",
                                    "allother"))


splitoutold <- read.csv("pyindex/2022/allagg_fall_500_lennosst_ALLsplit_biascorrect/Index.csv")
splitoutnew <- read.csv("pyindex/allagg_fall_500_lennosst_ALLsplit_biascorrect/Index.csv")

splitoutput2022 <- splitoutold %>%
  dplyr::mutate(Season = "fall",
                modname = "Survey10Prey") |>
  dplyr::left_join(stratlook) %>%
  dplyr::filter(Region %in% c("GOM", "GB", "MAB","MABGBstate", "bfin")) %>%
  dplyr::mutate(Type = ifelse(Region %in% c("GOM", "GB", "MAB"), "Ecoregion", "Bluefish"),
                Region = case_when(Region == "MABGBstate" ~ "StateWaters",
                                   Region == "bfin" ~ "SurveyBluefish",
                                   TRUE ~ Region))

splitoutput2023 <- splitoutnew %>%
  dplyr::mutate(Season = "fall",
                modname = "Combined25Prey") |>
  dplyr::left_join(stratlook) %>%
  dplyr::filter(Region %in% c("GOM", "GB", "MAB","MABGBstate", "bfin")) %>%
  dplyr::mutate(Type = ifelse(Region %in% c("GOM", "GB", "MAB"), "Ecoregion", "Bluefish"),
                Region = case_when(Region == "MABGBstate" ~ "StateWaters",
                                   Region == "bfin" ~ "SurveyBluefish",
                                   TRUE ~ Region))

foragemax <- max(splitoutput2023$Estimate) #scale to new prey list

splitoutput <- splitoutput2022 |> bind_rows(splitoutput2023)

foragetsmean <- splitoutput %>%
  dplyr::group_by(modname, Region) %>%
  dplyr::mutate(fmean = mean(Estimate)) 

ggplot(foragetsmean |> filter(Season =="fall"), aes(x=Time, y=((Estimate-fmean)/fmean), colour=modname)) +
  #geom_errorbar(aes(ymin=(Estimate+Std..Error.for.Estimate - fmean)/fmean, 
  #                  ymax=(Estimate-Std..Error.for.Estimate - fmean)/fmean))+
  geom_point(aes(shape=modname))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  #scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Relative forage biomass scaled to time series mean")  +
  ggtitle("Fall models: trend comparison") +
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0),
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))

```

```{r preycut2, fig.cap="Scale comparison between fall forage indices using different bluefish prey cut-offs: 10 or more observations by survey and 25 observations combined across surveys (current method).", fig.height=8}

ggplot(splitoutput |> filter(Season =="fall"), aes(x=Time, y=Estimate, colour=modname)) +
  #geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
  geom_ribbon(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate, fill=modname), linetype = 0, alpha = 0.15)+
  geom_point(aes(shape=modname))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Relative forage biomass scaled to maximum")  +
  ggtitle("Fall models: scale comparison")+
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0),
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))


# fig <- ggplot(splitoutput, aes(x=Time, y=Estimate, colour=Region)) +
#   geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
#   geom_point()+
#   geom_line()+
#   facet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
#   scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
#   ylab("Relative forage biomass scaled to maximum") 
# 
# 
# 
# EPU <- ggplot(splitoutput %>% dplyr::filter(Region %in% c("GOM", "GB", "MAB")), 
#               aes(x=Time, y=Estimate, colour=Region)) +
#   geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
#   geom_point()+
#   geom_line()+
#   geom_line(data=splitoutput2022 %>% dplyr::filter(Region %in% c("GOM", "GB", "MAB")),
#             aes(x=Time, y=Estimate, colour=Region), linetype = "dashed") +
#   facet_wrap(~fct_relevel(Region, "GOM", "GB", "MAB"), scales = "free_y", ncol = 1) +
#   scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
#   ylab("Relative forage biomass scaled to maximum") 
# 
# BF <- ggplot(splitoutput %>% dplyr::filter(Region %in% c("StateWaters", "SurveyBluefish")), 
#              aes(x=Time, y=Estimate, colour=Region)) +
#   geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
#   geom_point()+
#   geom_line()+
#   geom_line(data=splitoutput2022 %>% dplyr::filter(Region %in% c("StateWaters", "SurveyBluefish")),
#             aes(x=Time, y=Estimate, colour=Region), linetype = "dashed") +
#   facet_wrap(~fct_relevel(Region, "StateWaters", "SurveyBluefish"), scales = "free_y", ncol = 1) +
#   scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
#   ylab("Relative forage biomass scaled to maximum")


```


## Exclude assessed prey (Atlantic herring, Atlantic mackerel, silver hake) {-}

While the index is not sensitive to slight changes in the cut-off for including rarely observed prey, we would expect sensitivity to inclusion or exclusion of major prey. For a different application evaluating only unmanaged forage species, three assessed prey were excluded from the dataset: Atlantic herring, Atlantic mackerel, and silver hake. In NEFSC stomach data, unidentified Clupeids are also likely to be Atlantic herring, so this prey group was excluded as well. Using this dataset, we can evaluate the sensitivity of the index to major contributors to the prey list. We hypothesized that excluding these three assessed prey might make the most difference at the larger spatial scales including offshore regions (EPUs, Gulf of Maine, Georges Bank, and Mid-Atlantic). In particular, the dominance of silver hake in stomach content inputs from the Gulf of Maine and of silver hake and herring in stomach content inputs on Georges Bank noted in the manuscript suggests that excluding these species in those regions should affect the index.

The NEAMAP dataset had few of these three species observed in predator stomach contents sampled nearshore in the Mid-Atlantic bight, therefore we hypothesized that the index might be less sensitive to inclusion or exclusion at the scale of the bluefish assessment indices in state waters and bluefish survey strata.

As hypothesized, there are differences in both scale and trend for prey indices estimated by VAST for the Gulf of Maine and Georges Bank when excluding assessed prey species (Figs. Sd\@ref(fig:noassessedprey1)-Sd\@ref(fig:noassessedprey2)). There is a difference in scale for the Mid Atlantic, but index trends correspond better in that region.

Estimated indices at the bluefish assessment scale (state waters and bluefish inshore survey strata) have remarkably similar trends, although leaving out the three assessed prey species reduces the scale of the index.


```{r noassessedprey1, fig.cap="Trend comparison between fall forage indices using different prey groups: exclude assessed prey (Atlantic herring, silver hake and Atlantic mackerel) vs. all bluefish prey (current method).", fig.height=8}
splitoutunmod <- read.csv("pyindex/unmod/allagg_fall_500_lennosst_ALLsplit_biascorrect/Index.csv")

splitoutputunmod <- splitoutunmod %>%
  dplyr::mutate(Season = "fall",
                modname = "noHerringHakeMack") |>
  dplyr::left_join(stratlook) %>%
  dplyr::filter(Region %in% c("GOM", "GB", "MAB","MABGBstate", "bfin")) %>%
  dplyr::mutate(Type = ifelse(Region %in% c("GOM", "GB", "MAB"), "Ecoregion", "Bluefish"),
                Region = case_when(Region == "MABGBstate" ~ "StateWaters",
                                   Region == "bfin" ~ "SurveyBluefish",
                                   TRUE ~ Region))
splitoutput2023 <- splitoutput2023 %>%
  dplyr::mutate(modname = "allPrey")

splitoutput <- splitoutputunmod |> bind_rows(splitoutput2023)

foragetsmean <- splitoutput %>%
  dplyr::group_by(modname, Region) %>%
  dplyr::mutate(fmean = mean(Estimate)) 

ggplot(foragetsmean |> filter(Season =="fall"), aes(x=Time, y=((Estimate-fmean)/fmean), colour=modname)) +
  #geom_errorbar(aes(ymin=(Estimate+Std..Error.for.Estimate - fmean)/fmean, 
  #                  ymax=(Estimate-Std..Error.for.Estimate - fmean)/fmean))+
  geom_point(aes(shape=modname))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  #scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Relative forage biomass scaled to time series mean")  +
  ggtitle("Fall models: trend comparison") +
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0),
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))

```

```{r noassessedprey2, fig.cap="Scale comparison between fall forage indices using different prey groups: exclude assessed prey (Atlantic herring, silver hake and Atlantic mackerel) vs. all bluefish prey (current method).", fig.height=8}

ggplot(splitoutput |> filter(Season =="fall"), aes(x=Time, y=Estimate, colour=modname)) +
  #geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
  geom_ribbon(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate, fill=modname), linetype = 0, alpha = 0.15)+
  geom_point(aes(shape=modname))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Relative forage biomass scaled to maximum")  +
  ggtitle("Fall models: scale comparison")+
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0),
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))
# 
# 
# EPU <- ggplot(splitoutput %>% dplyr::filter(Region %in% c("GOM", "GB", "MAB")), 
#               aes(x=Time, y=Estimate, colour=Region)) +
#   geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
#   geom_point()+
#   geom_line()+
#   geom_line(data=splitoutput2022 %>% dplyr::filter(Region %in% c("GOM", "GB", "MAB")),
#             aes(x=Time, y=Estimate, colour=Region), linetype = "dashed") +
#   geom_line(data=splitoutputunmod %>% dplyr::filter(Region %in% c("GOM", "GB", "MAB")),
#             aes(x=Time, y=Estimate, colour=Region), linetype = "dotted") +
#   facet_wrap(~fct_relevel(Region, "GOM", "GB", "MAB"), scales = "free_y", ncol = 1) +
#   scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
#   ylab("Relative forage biomass scaled to maximum") 
# 
# BF <- ggplot(splitoutput %>% dplyr::filter(Region %in% c("StateWaters", "SurveyBluefish")), 
#              aes(x=Time, y=Estimate, colour=Region)) +
#   geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
#   geom_point()+
#   geom_line()+
#   geom_line(data=splitoutput2022 %>% dplyr::filter(Region %in% c("StateWaters", "SurveyBluefish")),
#             aes(x=Time, y=Estimate, colour=Region), linetype = "dashed") +
#   geom_line(data=splitoutputunmod %>% dplyr::filter(Region %in% c("StateWaters", "SurveyBluefish")),
#             aes(x=Time, y=Estimate, colour=Region), linetype = "dotted") +
#   facet_wrap(~fct_relevel(Region, "StateWaters", "SurveyBluefish"), scales = "free_y", ncol = 1) +
#   scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
#   ylab("Relative forage biomass scaled to maximum")
# 
# EPU
# 
# BF

```


# Forage index sensitivity to included predator species {-}

Preliminary analysis demonstrated that the first step of model selection was robust to both the inclusion/exclusion of NEAMAP data and the use of a different piscivore list based on an older diet similarity analysis in @garrison_dietary_2000. The full model with anisotrophy, spatial random effects and spatio-temporal random effects was always selected (see [NEFSC old predator list](https://sgaichas.github.io/bluefishdiet/VASTmodsel.html) and [NEFSC/NEAMAP current predator list](https://sgaichas.github.io/bluefishdiet/VASTmodsel_updatedPreds.html)).

Based on previous work [@ng_predator_2021], we expect that individual predator stomach contents only partially reflect available forage. This was the reason for including multiple predators to construct the forage index, and selecting the subset of predator size classes with the most similar diet to our target predator, bluefish. Therefore, we do expect some index sensitivity to which predators are included. 

Similar to the prey inclusion evaluation above, we expect the index to be less sensitive to the inclusion of predators with relatively low sampling or sampling limited to a small proportion of years or areas relative to predators with relatively high sampling levels overall that have been sampled consistently through time and across space.

The number of stomachs currently included in the NEFSC dataset varies considerably by predator/size class. For example, while there are many Atlantic cod stomachs in the NEFSC database, we used only stomachs from extra large cod, the size class with the most similar diet to all sizes of bluefish. Based on our assumptions above, we would expect the index to be less sensitive to the removal of cusk, fourspot flounder, or one of the squid species that contribute less than 1% of stomachs to the total than it would be to the removal of spiny dogfish or silver hake which each contribute nearly 25% of stomachs to the total.

```{r}
# object is called `allfh`
load(here("fhdat/allfh.Rdata"))

#object is called allfh21
load(here("fhdat/allfh21.Rdata"))

allfh <- allfh %>%
  dplyr::bind_rows(allfh21)

dietoverlap <- read_csv(here("fhdat/tgmat.2022-02-15.csv"))

library(dendextend)

d_dietoverlap <- dist(dietoverlap)

guilds <- hclust(d_dietoverlap)

#plot(guilds)

dend <- as.dendrogram(guilds)

dend <- rotate(dend, 1:136)

dend <- color_branches(dend, k=6) # Brian uses 6 categories

labels(dend) <- paste(as.character(names(dietoverlap[-1]))[order.dendrogram(dend)],
                           "(",labels(dend),")", 
                           sep = "")

dend <- hang.dendrogram(dend,hang_height=0.1)

pisccomplete <- partition_leaves(dend)[[
  which_node(dend, c("Bluefish..S(37)", "Bluefish..M(36)", "Bluefish..L(35)"))
]]

pisccompleten <- str_remove(pisccomplete, "\\(.*\\)")

dietoverlappisc <- dietoverlap |> 
  rename(pred = "...1") |>
  filter(pred %in% pisccompleten)

# #stats on overlap with bluefish
# # all preds
# mean(dietoverlap$Bluefish..L)
# mean(dietoverlap$Bluefish..M)
# mean(dietoverlap$Bluefish..S)
# 
# 
# # just piscivore preds
# mean(dietoverlappisc$Bluefish..L)
# mean(dietoverlappisc$Bluefish..M)
# mean(dietoverlappisc$Bluefish..S)
# 
# summarise_all(dietoverlappisc|>select(-pred), mean) #list(min, mean, max))
# 

pisccompletedf <- data.frame("COMNAME" = toupper(str_remove(pisccomplete, "\\..*")),
                              "SizeCat" = str_remove(str_extract(pisccomplete, "\\..*[:upper:]+"), "\\.."),
                              "feedguild" = "pisccomplete")
 
 fh.nefsc.pisc.pisccomplete <- allfh %>%
  #filter(pynam != "EMPTY") %>%
  left_join(pisccompletedf, by = c("pdcomnam" = "COMNAME",
                               "sizecat" = "SizeCat")) %>%
  filter(!is.na(feedguild))
 
NEAMAPdat <- read.csv(here("fhdat/NEAMAP_Bluefish Prey Index_Pred Stomach Counts.csv"))

NEAMAPstomn <- NEAMAPdat |>
  pivot_longer(-year, names_to = "comnam", values_to = "pdcount") |>
  mutate(pdcomnam = str_replace(str_to_upper(comnam), "\\.", " "))

 
stomcount <- fh.nefsc.pisc.pisccomplete |> 
  filter(year > 1984) |>
  group_by(cruise6, pdcomnam, station, pdid) |> 
  summarise(pycount = n()) |>
  group_by(pdcomnam) |>
  summarise(pdcount = n()) |>
  left_join(NEAMAPstomn |> group_by(pdcomnam) |> summarize(neamapcount = sum(pdcount))) |>
  rowwise() |>
  mutate('N stomachs' = sum(c_across(pdcount:neamapcount), na.rm = T)) |>
  select(Predator = pdcomnam, 'N stomachs')

flextable::flextable(stomcount) %>%
  flextable::set_caption("Total number of stomachs for each predator included in analysis, 1985-2021") %>%
  flextable::fontsize(size = 9, part = "all") %>%
  flextable::autofit()
```

The number of stomachs by year for each predator show general stability that seem unlikely to explain fluctuations in the index. Tables Sd\@ref(tab:nstom1) - Sd\@ref(tab:nstom11) show number of stomachs by predator included from both the NEFSC and NEAMAP datasets. While there was more sampling overall during the 1990s relative to other portions of the time series, variation in the full Northeast US annual forage index does not correspond to higher or lower numbers of predator stomachs included as input data (Fig. Sd\@ref(fig:nstom2)).

```{r}
block_section(op_section)
```


```{r nstom1, ft.arraystretch = 1}
stomcountyr1 <- fh.nefsc.pisc.pisccomplete |>
  mutate(year = as.numeric(year)) |>
  filter(year %in% c(1984:2001)) |>
  group_by(cruise6, year, pdcomnam, station, pdid) |> 
  summarise(pycount = n()) |>
  group_by(pdcomnam, year) |>
  summarise(pdcount = n()) |>
  #bind_rows(NEAMAPstomn |> select(-comnam)) |>
  group_by(pdcomnam, year) |>
  summarize(pdcount = sum(pdcount, na.rm = T)) |>
  arrange(year, pdcomnam) |>
  pivot_wider(names_from = year, values_from = pdcount)

flextable::flextable(stomcountyr1) %>%
  flextable::set_header_labels(pdcomnam = "predator") |>
  flextable::set_caption("Number of stomachs for each predator and year included in analysis, 1985-2001") %>%
  #flextable::rotate(j = 2:ncol_keys(.), align = "bottom", rotation = "btlr", part = "header") |>
  #flextable::height(height = 1.2, part = "header") |>
  flextable::padding(padding = 0, part = "all") %>%
  flextable::fontsize(size = 7, part = "all") %>%
  #flextable::width(width = dim_pretty((.), part = "all", unit = "in", hspans = "none")) 
  flextable::set_table_properties(layout = "autofit", width = 1) |>
  flextable::theme_zebra()

```

```{r nstom11, ft.arraystretch = 1}
stomcountyr2 <- fh.nefsc.pisc.pisccomplete |>
  mutate(year = as.numeric(year)) |>
  filter(year > 2001) |>
  group_by(cruise6, year, pdcomnam, station, pdid) |> 
  summarise(pycount = n()) |>
  group_by(pdcomnam, year) |>
  summarise(pdcount = n()) |>
  bind_rows(NEAMAPstomn |> select(-comnam)) |>
  group_by(pdcomnam, year) |>
  summarize(pdcount = sum(pdcount, na.rm = T)) |>
  arrange(year, pdcomnam) |>
  pivot_wider(names_from = year, values_from = pdcount)

flextable::flextable(stomcountyr2) %>%
  flextable::set_header_labels(pdcomnam = "predator") |>
  flextable::set_caption("Number of stomachs for each predator and year included in analysis, 2002-2021") %>%
  #flextable::rotate(j = 2:ncol_keys(.), align = "bottom", rotation = "btlr", part = "header") |>
  #flextable::height(height = 1.2, part = "header") |>
  flextable::padding(padding = 0, part = "all") %>%
  flextable::fontsize(size = 7, part = "all") %>%
  #flextable::width(width = dim_pretty((.), part = "all", unit = "in", hspans = "none")) 
  flextable::set_table_properties(layout = "autofit", width = 1) |>
  flextable::theme_zebra()

```


```{r}
block_section(close_section)
```


```{r nstom2, fig.cap="Full Northeast US shelf annual forage index compared with total number of stomachs used as input."}

stomcountyr <- fh.nefsc.pisc.pisccomplete |>
  filter(year > 1984) |>
  group_by(cruise6, year, pdcomnam, station, pdid) |> 
  summarise(pycount = n()) |>
  group_by(pdcomnam, year) |>
  summarise(pdcount = n()) |>
  bind_rows(NEAMAPstomn |> select(-comnam)) |>
  group_by(pdcomnam, year) |>
  summarize(pdcount = sum(pdcount, na.rm = T)) |>
  arrange(year, pdcomnam) |>
  pivot_wider(names_from = year, values_from = pdcount)


outann <- read.csv("pyindex/allagg_annual_500_lennosst_ALLsplit_biascorrect/Index.csv")

splitoutann <-outann %>%
  dplyr::left_join(stratlook) %>%
  dplyr::filter(Region %in% c("AllEPU", "GOM", "GB", "MAB","MABGBstate", "bfin")) %>%
  dplyr::mutate(Type = ifelse(Region %in% c("AllEPU","GOM", "GB", "MAB"), "Ecoregion", "Bluefish"),
                Region = case_when(Region == "MABGBstate" ~ "StateWaters",
                                   Region == "bfin" ~ "SurveyBluefish",
                                   TRUE ~ Region))

foragemaxann <- max(splitoutann$Estimate, na.rm = T)

AllEPU <- ggplot(splitoutann %>% dplyr::filter(Region %in% c("AllEPU")), 
              aes(x=Time, y=Estimate, colour=Region)) +
  geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
  geom_point()+
  geom_line() +
  scale_y_continuous(labels=function(x)round(x/foragemaxann, digits = 2)) +
  ylab("Relative forage biomass scaled to maximum") +
  ggtitle("Annual forage index, full Northeast US shelf (GOM + GB + MAB)")

NEAMAPtot <- read.csv(here("fhdat/NEAMAP_Mean stomach weights_Bluefish Prey_wWQ2023.csv"))

NEAMAPtotstomn <- NEAMAPtot |>
  select(year, season, nstomtot) |>
  group_by(year) |>
  summarise(pdcount = sum(nstomtot)) |>
  mutate(pdcomnam = "NEAMAP")

stackbar <- stomcountyr |>
  pivot_longer(-pdcomnam, names_to = "year", values_to = "pdcount") |>
  mutate(year = as.integer(year)) |>
  # fh.nefsc.pisc.pisccomplete |> 
  # filter(year > 1984) |>
  # group_by(cruise6, year, pdcomnam, station, pdid) |> 
  # summarise(pycount = n()) |>
  # group_by(pdcomnam, year) |>
  # summarise(pdcount = n()) |>
  # ungroup() |>
  # arrange(year, pdcomnam) |>
  # bind_rows(NEAMAPstomn) |>
  ggplot(aes(x=year, y=pdcount, fill=pdcomnam)) +
  geom_bar(position = "stack", stat = "identity") +
  theme(legend.title = element_blank(), 
        legend.text = element_text(size = 6),
        legend.key.size = unit(0.3, "cm")) +
  guides(fill = guide_legend(ncol = 1)) +
  viridis::scale_fill_viridis(option="turbo", discrete=TRUE) +
  ggtitle("Number of stomachs included in dataset by year")


stomtotyr <- stomcountyr |>
  pivot_longer(-pdcomnam, names_to = "year", values_to = "pdcount") |>
  mutate(year = as.integer(year)) |> 
  group_by(year) |> 
  summarise(totsum = sum(pdcount, na.rm=T)) |> 
  as.data.frame()

#stackbar

#AllEPU 

splitoutann %>% dplyr::filter(Region %in% c("AllEPU")) |>
  ggplot(aes(x=Time, y=Estimate)) +
  geom_line() +
  geom_line(aes(x=Time, y=mean(Estimate)), linetype=2) +
  #geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate)) +
  geom_ribbon(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate), alpha = 0.2) + 
  geom_line(data=stomtotyr, aes(x=year, y=totsum), color="blue") +
  geom_line(data=stomtotyr, aes(x=year, y=mean(totsum)), color="blue", linetype=2) +
  ylab("Number of stomachs/estimated annual forage")  +
  ggtitle("Total stomachs compared with full shelf annual index")+
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0),
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))
            
```


Here we illustrate a range of sensitivities from removing four predators individually and examining the impact on both model convergence and impact to the index.

```{r getmodinfoindex}

# from each output folder in pyindex, 
outdir <- here("pyindex")
moddirs <- list.dirs(outdir)
# only get pred sensitivities and original full models
moddirssens <- grep("predsens", moddirs, value = TRUE)
fullfall <- grep("pyindex/allagg_fall_", moddirs, value = TRUE)
fullspring <- grep("pyindex/allagg_spring_", moddirs, value = TRUE)
moddirs <- c(moddirssens, fullfall, fullspring)
# keep folder name
modnames <- list.dirs(outdir, full.names = FALSE)


# function to apply extracting info
getmodinfo <- function(d.name){
  # read settings
  modpath <- stringr::str_split(d.name, "/", simplify = TRUE)
  modname <- modpath[length(modpath)]
  
  settings <- read.table(file.path(d.name, "settings.txt"), comment.char = "",
    fill = TRUE, header = FALSE)
  
  n_x <- as.numeric(as.character(settings[(which(settings[,1]=="$n_x")+1),2]))
  grid_size_km <- as.numeric(as.character(settings[(which(settings[,1]=="$grid_size_km")+1),2]))
  max_cells <- as.numeric(as.character(settings[(which(settings[,1]=="$max_cells")+1),2]))
  use_anisotropy <- as.character(settings[(which(settings[,1]=="$use_anisotropy")+1),2])
  fine_scale <- as.character(settings[(which(settings[,1]=="$fine_scale")+1),2])
  bias.correct <- as.character(settings[(which(settings[,1]=="$bias.correct")+1),2])
  
  #FieldConfig
  if(settings[(which(settings[,1]=="$FieldConfig")+1),1]=="Component_1"){
    omega1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+2),2])
    omega2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+3),1])
    epsilon1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+4),2])
    epsilon2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+5),1])
    beta1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+6),2])
    beta2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+7),1])
  }
  
  if(settings[(which(settings[,1]=="$FieldConfig")+1),1]=="Omega1"){
    omega1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+3),1])
    omega2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+4),1])
    epsilon1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+3),2])
    epsilon2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+4),2])
    beta1 <- "IID"
    beta2 <- "IID"
  }
  
  
  #RhoConfig
  rho_beta1 <- as.numeric(as.character(settings[(which(settings[,1]=="$RhoConfig")+3),1]))
  rho_beta2 <- as.numeric(as.character(settings[(which(settings[,1]=="$RhoConfig")+3),2]))
  rho_epsilon1 <- as.numeric(as.character(settings[(which(settings[,1]=="$RhoConfig")+4),1]))
  rho_epsilon2 <- as.numeric(as.character(settings[(which(settings[,1]=="$RhoConfig")+4),2]))
  
  # read parameter estimates, object is called parameter_Estimates
  load(file.path(d.name, "parameter_estimates.RData"))
  
  AIC <- parameter_estimates$AIC[1]  
  converged <- parameter_estimates$Convergence_check[1]
  fixedcoeff <- unname(parameter_estimates$number_of_coefficients[2])
  randomcoeff <- unname(parameter_estimates$number_of_coefficients[3])
  
  #index <- read.csv(file.path(d.name, "Index.csv"))
  
  # return model attributes as a dataframe
  out <- data.frame(modname = modname,
                    n_x = n_x,
                    grid_size_km = grid_size_km,
                    max_cells = max_cells,
                    use_anisotropy = use_anisotropy,
                    fine_scale =  fine_scale,
                    bias.correct = bias.correct,
                    omega1 = omega1,
                    omega2 = omega2,
                    epsilon1 = epsilon1,
                    epsilon2 = epsilon2,
                    beta1 = beta1,
                    beta2 = beta2,
                    rho_epsilon1 = rho_epsilon1,
                    rho_epsilon2 = rho_epsilon2,
                    rho_beta1 = rho_beta1,
                    rho_beta2 = rho_beta2,
                    AIC = AIC,
                    converged = converged,
                    fixedcoeff = fixedcoeff,
                    randomcoeff = randomcoeff#,
                    #index = index
  )
  	
	return(out)

}

# getfitdat <- function(d.name){
#   # read settings
#   modpath <- stringr::str_split(d.name, "/", simplify = TRUE)
#   modname <- modpath[length(modpath)]
#   
#   readRDS(file.path(d.name, "fit.rds"))
#   nstations <- fit$data_list$n_i # this is grid size not input N stations
#   out <- data.frame(modname = modname,
#                     nstations = nstations
#   )
#   
#   return(out)
# }


#Get number of stations in input data for each dataset

predout <- c("", "_no4spot", "_nololigo", "_nowhake", "_nospdog")

modnstations <- data.frame()

for(dat in predout){
  
  # this dataset created in SSTmethods.Rmd
  
  fname <- paste0("fhdat/bluepyagg_stn_all_OISST", dat, ".rds")
  #print(fname)
  
  #bluepyagg_stn <- readRDS(here::here("fhdat/bluepyagg_stn_all_OISST.rds"))
  bluepyagg_stn <- readRDS(here::here(fname))
  
  # make SST column that uses surftemp unless missing or 0
  # there are 3 surftemp 0 values in the dataset, all with oisst > 15
  bluepyagg_stn <- bluepyagg_stn %>%
    dplyr::mutate(sstfill = ifelse((is.na(surftemp)|surftemp==0), oisst, surftemp))
  
  #bluepyagg_stn <- bluepyagg_pred_stn
  
  # filter to assessment years at Tony's suggestion
  
  # code Vessel as AL=0, HB=1, NEAMAP=2
  
  bluepyagg_stn_fall <- bluepyagg_stn %>%
    #ungroup() %>%
    filter(season_ng == "FALL",
           year > 1984) %>%
    mutate(AreaSwept_km2 = 1, #Elizabeth's code
           #declon = -declon already done before neamap merge
           Vessel = as.numeric(as.factor(vessel))-1
    ) %>% 
    dplyr::select(Catch_g = meanbluepywt, #use bluepywt for individuals input in example
                  Year = year,
                  Vessel,
                  AreaSwept_km2,
                  Lat = declat,
                  Lon = declon,
                  meanpisclen,
                  npiscsp,
                  #bottemp, #this leaves out many stations for NEFSC
                  #surftemp, #this leaves out many stations for NEFSC
                  oisst,
                  sstfill
    ) %>%
    na.omit() %>%
    as.data.frame()
  
  
  bluepyagg_stn_spring <- bluepyagg_stn %>%
    filter(season_ng == "SPRING",
           year > 1984) %>%
    mutate(AreaSwept_km2 =1, #Elizabeth's code
           #declon = -declon already done before neamap merge
           Vessel = as.numeric(as.factor(vessel))-1
    ) %>% 
    dplyr::select(Catch_g = meanbluepywt,
                  Year = year,
                  Vessel,
                  AreaSwept_km2,
                  Lat = declat,
                  Lon = declon,
                  meanpisclen,
                  npiscsp,
                  #bottemp, #this leaves out many stations for NEFSC
                  #surftemp, #this leaves out many stations for NEFSC
                  oisst,
                  sstfill
    ) %>%
    na.omit() %>%
    as.data.frame()
  
  nstations <- data.frame(dat = map_chr(str_split(dat, "_"), 2, .default = NA_character_),
                          season = c("fall", "spring"),
                          nstations = c(dim(bluepyagg_stn_fall)[1],
                                        dim(bluepyagg_stn_spring)[1])
  )
  
  modnstations <- bind_rows(modnstations, nstations)
}


# function to apply extracting info
getmodindex <- function(d.name){
  # read settings
  modpath <- stringr::str_split(d.name, "/", simplify = TRUE)
  modname <- modpath[length(modpath)]
  
  index <- read.csv(file.path(d.name, "Index.csv"))
  # return model indices as a dataframe
  out <- data.frame(modname = modname,
                    index
  )
  
  return(out)
}

modcompare <- purrr::map_dfr(moddirs, getmodinfo)

modcompareindex <- purrr::map_dfr(moddirs, getmodindex)

```

All of the models converged with alternative predator inputs. We report the total number of stations included for each model (nstations) and the number of stations lost relative to the full dataset with all predators (stationloss; Table Sd\@ref(tab:altpreds)). The input datasets differ, so comparing AIC across these models is not appropriate, but we report AIC by model for information.

```{r altpreds}

modsense <- modcompare %>%
  dplyr::mutate(season = case_when(str_detect(modname, "_fall_") ~ "Fall",
                            str_detect(modname, "spring") ~ "Spring",
                            str_detect(modname, "_all_") ~ "Annual",
                            TRUE ~ as.character(NA))) %>%
  dplyr::mutate(dat = modname |> map(str_split, pattern = "_") |> map_chr(c(1,7), .default = NA_character_)) %>%
  dplyr::mutate(converged2 = case_when(str_detect(converged, "no evidence") ~ "likely",
                                str_detect(converged, "is likely not") ~ "unlikely",
                                TRUE ~ as.character(NA))) %>%
  dplyr::left_join(modnstations |> mutate(season = str_to_sentence(season))) |>
  dplyr::group_by(season) %>%
  dplyr::mutate(stationloss = max(nstations) - nstations,
                dat = if_else(is.na(dat), "allpreds", dat)) |>
  #dplyr::mutate(deltaAIC = AIC-min(AIC)) %>%
  dplyr::select(dat, season, #deltaAIC, 
                fixedcoeff,
         randomcoeff, 
         #use_anisotropy, 
         #omega1, omega2, epsilon1, epsilon2, 
         #beta1, beta2, 
         AIC, converged2, 
         nstations,
         stationloss) %>%
  dplyr::arrange(season, AIC)

# DT::datatable(modselect2, rownames = FALSE, 
#               options= list(pageLength = 50, scrollX = TRUE),
#               caption = "Comparison of delta AIC values for alternative vessel effects and catchability covariates. See text for model descriptions.")

flextable::flextable(modsense) %>%
  flextable::set_caption("Comparison of model convergence and sample size for predator sensitivity runs. See text for model descriptions.") %>%
  flextable::fontsize(size = 9, part = "all") %>%
  flextable::autofit()
```

## Impact of removing predators contributing fewer stomachs {-}

Removing fourspot flounder (demersal predator, small sample) and longfin squid (pelagic predator, small sample) resulted in the same or a similar number of survey stations included in analysis, with almost no change to the fall and spring indices trends (Figs. Sd\@ref(fig:comptrend1)-Sd\@ref(fig:comptrend2)) or scale (Figs. Sd\@ref(fig:compscale1)-Sd\@ref(fig:compscale2)). Squid stomachs were sampled only prior to the year 2000, but no change for the early portion of the index is apparent when leaving this predator out.

```{r comptrend1, fig.cap="Trend comparison between fall forage indices using different predator groups: exclude predators with low sampling (fourspot flounder, longfin squid) vs. full predator list (current method).", fig.height=8}
splitoutput <- modcompareindex %>%
  dplyr::mutate(Season = modname |> map(str_split, pattern = "_") |> map_chr(c(1,2))) %>%
  dplyr::mutate(Data = modname |> map(str_split, pattern = "_") |> map_chr(c(1,7), .default = "allpreds")) %>%
  dplyr::left_join(stratlook) %>%
  dplyr::filter(Region %in% c("GOM", "GB", "MAB","MABGBstate", "bfin")) %>%
  dplyr::mutate(Type = ifelse(Region %in% c("GOM", "GB", "MAB"), "Ecoregion", "Bluefish"),
                Region = case_when(Region == "MABGBstate" ~ "StateWaters",
                                   Region == "bfin" ~ "SurveyBluefish",
                                   TRUE ~ Region))

foragemax <- max(splitoutput$Estimate)


foragetsmean <- splitoutput %>%
  dplyr::group_by(modname, Region) %>%
  dplyr::mutate(fmean = mean(Estimate)) 


ggplot(foragetsmean |> filter(Season =="fall",
                             Data %in% c("allpreds",
                                         "no4spot",
                                         "nololigo")), 
       aes(x=Time, y=((Estimate-fmean)/fmean), colour=Data)) +
  #geom_errorbar(aes(ymin=(Estimate+Std..Error.for.Estimate - fmean)/fmean, 
  #                  ymax=(Estimate-Std..Error.for.Estimate - fmean)/fmean))+
  geom_point(aes(shape=Data))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  #scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Relative forage biomass scaled to time series mean")  +
  ggtitle("Fall models: trend comparison") +
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0)
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))

```

```{r comptrend2, fig.cap="Trend comparison between spring forage indices using different predator groups: exclude predators with low sampling (fourspot flounder, longfin squid) vs. full predator list (current method).", fig.height=8}

ggplot(foragetsmean |> filter(Season =="spring",
                             Data %in% c("allpreds",
                                         "no4spot",
                                         "nololigo")), 
       aes(x=Time, y=((Estimate-fmean)/fmean), colour=Data)) +
  #geom_errorbar(aes(ymin=(Estimate+Std..Error.for.Estimate - fmean)/fmean,
  #                  ymax=(Estimate-Std..Error.for.Estimate - fmean)/fmean))+
  geom_point(aes(shape=Data))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  #scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Relative forage biomass scaled to time series mean")  +
  ggtitle("Spring models: trend comparison")+
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0)
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))


```

```{r compscale1, fig.cap="Scale comparison between fall forage indices using different predator groups: exclude predators with low sampling (fourspot flounder, longfin squid) vs. full predator list (current method).", fig.height=8}


ggplot(splitoutput |> filter(Season =="fall",
                             Data %in% c("allpreds",
                                         "no4spot",
                                         "nololigo")), 
       aes(x=Time, y=Estimate, colour=Data)) +
  #geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
  geom_ribbon(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate, fill=Data), linetype = 0, alpha = 0.15)+
  geom_point(aes(shape=Data))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Relative forage biomass scaled to maximum")  +
  ggtitle("Fall models: scale comparison")+
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0)
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))

```

```{r compscale2, fig.cap="Scale comparison between spring forage indices using different predator groups: exclude predators with low sampling (fourspot flounder, longfin squid) vs. full predator list (current method).", fig.height=8}
  

ggplot(splitoutput |> filter(Season =="spring",
                             Data %in% c("allpreds",
                                         "no4spot",
                                         "nololigo")), 
       aes(x=Time, y=Estimate, colour=Data)) +
  #geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
  geom_ribbon(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate, fill=Data), linetype = 0, alpha = 0.15)+
  geom_point(aes(shape=Data))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Relative forage biomass scaled to maximum")  +
  ggtitle("Spring models: scale comparison")+
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0)
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))


```

## Impact of removing predators contributing many stomachs {-}

Removing a predator with more stomachs contributing to the input data, white hake, made small changes in the index, but well within the error bars of the original (Figs. Sd\@ref(fig:comptrend3)-Sd\@ref(fig:comptrend4), Sd\@ref(fig:compscale3)-Sd\@ref(fig:compscale4)). Changes in trend were most apparent in the Gulf of Maine.

Removing spiny dogfish, which contributes over 25% of stomachs in the analysis, has the largest impact on the index, as predicted. The primary impact is on index scale, which is lower when excluding spiny dogfish.

However, the forage index trend across multiple scales is remarkably robust to removing both white hake and spiny dogfish individually.


```{r comptrend3, fig.cap="Trend comparison between fall forage indices using different predator groups: exclude predators with high sampling (white hake, spiny dogfish) vs. full predator list (current method).", fig.height=8}


ggplot(foragetsmean |> filter(Season =="fall",
                             Data %in% c("allpreds",
                                         "nowhake",
                                         "nospdog")), 
       aes(x=Time, y=((Estimate-fmean)/fmean), colour=Data)) +
  #geom_errorbar(aes(ymin=(Estimate+Std..Error.for.Estimate - fmean)/fmean, 
  #                  ymax=(Estimate-Std..Error.for.Estimate - fmean)/fmean))+
  geom_point(aes(shape=Data))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Relative forage biomass scaled to time series mean")  +
  ggtitle("Fall models: trend comparison")+
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0)
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))

```


```{r comptrend4, fig.cap="Trend comparison between spring forage indices using different predator groups: exclude predators with high sampling (white hake, spiny dogfish) vs. full predator list (current method).", fig.height=8}
  

ggplot(foragetsmean |> filter(Season =="spring",
                             Data %in% c("allpreds",
                                         "nowhake",
                                         "nospdog")), 
       aes(x=Time, y=((Estimate-fmean)/fmean), colour=Data)) +
  #geom_errorbar(aes(ymin=(Estimate+Std..Error.for.Estimate - fmean)/fmean, 
  #                  ymax=(Estimate-Std..Error.for.Estimate - fmean)/fmean))+
  geom_point(aes(shape=Data))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Relative forage biomass scaled to time series mean")  +
  ggtitle("Spring models: trend comparison")+
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0)
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))


```


```{r compscale3, fig.cap="Scale comparison between fall forage indices using different predator groups: exclude predators with high sampling (white hake, spiny dogfish) vs. full predator list (current method).", fig.height=8}

ggplot(splitoutput |> filter(Season =="fall",
                             Data %in% c("allpreds",
                                         "nowhake",
                                         "nospdog")), 
       aes(x=Time, y=Estimate, colour=Data)) +
  #geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
  geom_ribbon(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate, fill=Data), linetype = 0, alpha = 0.15)+
  geom_point(aes(shape=Data))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Relative forage biomass scaled to maximum")  +
  ggtitle("Fall models: scale comparison")+
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0)
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))

```

```{r compscale4, fig.cap="Scale comparison between spring forage indices using different predator groups: exclude predators with high sampling (white hake, spiny dogfish) vs. full predator list (current method).", fig.height=8}

ggplot(splitoutput |> filter(Season =="spring",
                             Data %in% c("allpreds",
                                         "nowhake",
                                         "nospdog")), 
       aes(x=Time, y=Estimate, colour=Data)) +
  #geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
  geom_ribbon(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate, fill=Data), linetype = 0, alpha = 0.15)+
  geom_point(aes(shape=Data))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Relative forage biomass scaled to maximum")  +
  ggtitle("Spring models: scale comparison")+
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0)
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))


```


# VAST model fit and index sensitivity with depth as a covariate {-}

Previous work found predator size to be an important "catchability" covariate for stomach content based indices of Atlantic herring [@ng_predator_2021]. In the current work, the number of predator species was also found to be an important catchability covariate, along with sea surface temperature. Here we explore an additional potential covariate, depth, which might affect availability of pelagic forage species to demersal predators sampled by a bottom trawl.

Comparisons are presented for models with and without depth as a covariate. This sensitivity was run for the fall and spring models only, and include the "base" model (no covariates) and the two best fit models across all seasons (mean predator length, number of predators, and sea surface temperature as covariates, "lennumsst"; and mean predator length and number of predators as covariates "lennum") compared with depth only as a covariate, and depth combined with the two best fit models across all seasons: "deplennum" and "deplennumsst". 

There are 10 stations with depth information missing, which were dropped from this analysis. This leaves 26,590 stations combined over spring and fall.

```{r}

# from each output folder in pyindex, 
outdir <- here("pyindex_modsel2depth")
moddirs <- list.dirs(outdir) 
moddirs <- moddirs[-1]
# keep folder name
modnames <- list.dirs(outdir, full.names = FALSE)


modcompare <- purrr::map_dfr(moddirs, getmodinfo)
```


For model comparisons with different combinations of catchability covariates, "modname" in Table Sd\@ref(tab:modsel2) follows a similar pattern as in Supplement 2, with  all prey aggregated ("allagg" for all models), season ("spring" for models of months 1-6, "fall" for models of months 7-12), number of knots (500 for all models), and which covariates were included. The names for model options and associated VAST model settings are: 

Model covariates sensitivity options, FieldConfig default (all IID), with anisotropy:

1.  "_base" =          No vessel overdispersion or length/number covariates 
<!--
1.  "_len" =           Predator mean length covariate
1.  "_num" =           Number of predator species covariate
-->
1.  "_lennum" =        Predator mean length and number of predator species covariates
<!--
1.  "_sst" =           Combined in situ and OISST covariate
1.  "_lensst" =        Predator mean length and SST covariates
1.  "_numsst" =        Number of predator species and SST covariates
-->
1.  "_lennumsst" =     Predator mean length, number of predator species, and SST covariates
<!--
1.  "_eta10" =         Overdispersion (vessel effect) in first linear predictor (prey encounter)
1.  "_eta11" =         Overdispersion (vessel effect) in both linear predictors (prey encounter and weight)
-->
1.  "_depth" =         Survey station depth covariate
1.  "_deplennum" =     Survey station depth, predator mean length and number of predator species covariates
1.  "_deplennumsst" =  Survey station depth, predator mean length, number of predator species, and SST covariates


```{r modsel2}
# from each output folder in pyindex, 
outdir <- here("pyindex_modsel2depth")
moddirs <- list.dirs(outdir) 
moddirs <- moddirs[-1]
# keep folder name
modnames <- list.dirs(outdir, full.names = FALSE)

modcompare2 <- purrr::map_dfr(moddirs, getmodinfo)

modselect2 <- modcompare2 %>%
  dplyr::mutate(season = case_when(str_detect(modname, "_fall_") ~ "Fall",
                            str_detect(modname, "spring") ~ "Spring",
                            str_detect(modname, "_all_") ~ "Annual",
                            TRUE ~ as.character(NA))) %>%
  dplyr::mutate(converged2 = case_when(str_detect(converged, "no evidence") ~ "likely",
                                str_detect(converged, "is likely not") ~ "unlikely",
                                TRUE ~ as.character(NA))) %>%
  dplyr::group_by(season) %>%
  dplyr::mutate(deltaAIC = AIC-min(AIC)) %>%
  dplyr::select(modname, season, deltaAIC, fixedcoeff,
         randomcoeff, use_anisotropy, 
         omega1, omega2, epsilon1, epsilon2, 
         beta1, beta2, AIC, converged2) %>%
  dplyr::arrange(season, AIC)

# DT::datatable(modselect2, rownames = FALSE, 
#               options= list(pageLength = 50, scrollX = TRUE),
#               caption = "Comparison of delta AIC values for alternative vessel effects and catchability covariates. See text for model descriptions.")

flextable::flextable(modselect2%>%
                       dplyr::select(-c(use_anisotropy, 
         omega1, omega2, epsilon1, epsilon2, 
         beta1, beta2))) %>%
  flextable::set_caption("Comparison of delta AIC values for alternative vessel effects and catchability covariates. See text for model descriptions.") %>%
  flextable::fontsize(size = 9, part = "all") %>%
  flextable::autofit()


```


Model comparisons above for fall and spring reinforce initial model selection findings of the best model fit including mean predator length, number of predator species, and SST at a survey station as catchability covariates. However, models with depth included as a covariate had convergence issues, with maximum gradients several orders of magnitude higher than models without depth.


```{r indexcompare}
# 
# 
# modcompareindex <- purrr::map_dfr(moddirs, getmodindex)
# 
# splitoutput <- modcompareindex %>%
#   dplyr::mutate(Season = modname |> map(str_split, pattern = "_") |> map_chr(c(1,2))) %>%
#   dplyr::left_join(stratlook) %>%
#   dplyr::filter(Region %in% c("GOM", "GB", "MAB","MABGBstate", "bfin")) %>%
#   dplyr::mutate(Type = ifelse(Region %in% c("GOM", "GB", "MAB"), "Ecoregion", "Bluefish"),
#                 Region = case_when(Region == "MABGBstate" ~ "StateWaters",
#                                    Region == "bfin" ~ "SurveyBluefish",
#                                    TRUE ~ Region))
# 
# foragemax <- max(splitoutput$Estimate)
# 
# ggplot(splitoutput |> filter(Season =="fall"), aes(x=Time, y=Estimate, colour=modname)) +
#   geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
#   geom_point()+
#   geom_line()+
#   #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
#   facet_wrap(~Region, scales = "free_y") +
#   scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
#   ylab("Relative forage biomass scaled to maximum")  +
#   ggtitle("Fall models: scale comparison")+
#   theme(legend.position = c(1, 0),
#         legend.justification = c(1, 0),
#         legend.text = element_text(size=rel(0.5)),
#         legend.key.size = unit(0.5, 'lines'),
#         legend.title = element_text(size=rel(0.5)))
#   
# 
# # ggplot(splitoutput |> filter(Season =="spring"), aes(x=Time, y=Estimate, colour=modname)) +
# #   geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
# #   geom_point()+
# #   geom_line()+
# #   #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
# #   facet_wrap(~Region, scales = "free_y") +
# #   scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
# #   ylab("Relative forage biomass scaled to maximum")  +
# #   ggtitle("Spring models")
# 
# 
```


```{r}
# 
# foragetsmean <- splitoutput %>%
#   dplyr::group_by(modname, Region) %>%
#   dplyr::mutate(fmean = mean(Estimate)) 
# 
# 
# ggplot(foragetsmean |> filter(Season =="fall"), aes(x=Time, y=((Estimate-fmean)/fmean), colour=modname)) +
#   #geom_errorbar(aes(ymin=(Estimate+Std..Error.for.Estimate - fmean)/fmean, 
#   #                  ymax=(Estimate-Std..Error.for.Estimate - fmean)/fmean))+
#   geom_point()+
#   geom_line()+
#   #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
#   facet_wrap(~Region, scales = "free_y") +
#   #scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
#   ylab("Relative forage biomass scaled to time series mean")  +
#   ggtitle("Fall models: trend comparison") +
#   theme(legend.position = c(1, 0),
#         legend.justification = c(1, 0),
#         legend.text = element_text(size=rel(0.5)),
#         legend.key.size = unit(0.5, 'lines'),
#         legend.title = element_text(size=rel(0.5)))
# 
# # ggplot(foragetsmean |> filter(Season =="spring"), aes(x=Time, y=((Estimate-fmean)/fmean), colour=modname)) +
# #   #geom_errorbar(aes(ymin=(Estimate+Std..Error.for.Estimate - fmean)/fmean, 
# #   #                  ymax=(Estimate-Std..Error.for.Estimate - fmean)/fmean))+
# #   geom_point()+
# #   geom_line()+
# #   #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
# #   facet_wrap(~Region, scales = "free_y") +
# #   #scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
# #   ylab("Relative forage biomass scaled to time series mean")  +
# #   ggtitle("Spring models")
# 

```


# References {-}