Forage Index update for 2025 SOE

Data updated through 2023

1 Summary

I needed to compare the 1982-2022 submitted for the 2024 report with the 1982-2023 to be submitted for 2025 because a change to the VAST code caused my previous setup to break.

Also, its a full new input NEFSC food habits dataset, rather than tagging on years as was done in the past. NEAMAP just added one year from last year’s submission.

I can no longer produce the inshore <3 nmi and offshore strata which used a custom extrapolation grid. There may be a fix to estimate the bluefish assessment indices (see here), but the current run uses the built in grid and only estimates EPU level indices for the SOE.

2 Workflow

Operational updates to the forage index submitted to the State of the Ecosystem report are in the forageindex github repository: https://github.com/NOAA-EDAB/forageindex

Data input files are in the folder fhdat and were processed with the script VASTforage_ProcessInputDat.R in that folder: https://github.com/NOAA-EDAB/forageindex/blob/main/fhdat/VASTforage_ProcessInputDat.R

# Streamlined version of CreateVASTInputs.Rmd for operational updates to forage index indicators
# October 2023
#   This one is updating with 2022 NEFSC and NEAMAP data and OISST
#   To be used in the 2024 State of the Ecosystem report
# October 2024, added Brian's new input dataset, all else the same

library(tidyverse)
library(here)
library(dendextend)

# Load NEFSC stomach data received from Brian Smith

# New full dataset Septemer 2024 from https://drive.google.com/drive/folders/11aJkYoE1eVExKig-9rcMrgJ0V_GBIgkn
# downloaded datafile allfh.rmd.epu.Rdata and put it in my fhdat folder on this repo

# object is called `allfh`
load(here("fhdat/allfh.rmd.epu.Rdata"))

# Prior to Oct 2024, datasets were appended
# object is called `allfh`
# load(here("fhdat/allfh.Rdata"))
# 
# #object is called allfh21
# load(here("fhdat/allfh21.Rdata"))
# 
# #object is called allfh22
# load(here("fhdat/allfh22.Rdata"))
# 
# # bind all NEFSC stomach datasets
# allfh <- allfh %>%
#   dplyr::bind_rows(allfh21) |>
#   dplyr::bind_rows(allfh22)

###############################################################################
# 2024 note: diet overlap could be updated with new diet data, but 
# I DONT WANT THE PREDATOR LIST TO CHANGE SLIGHTLY EVERY YEAR WITHOUT TESTING
# Revisit after several years, do comparisons, then decide whether to update
#
# read predator similarity info to generate predator list
# Input NEFSC food habits overlap matrix:

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

# use dendextend functions to get list
d_dietoverlap <- dist(dietoverlap)

guilds <- hclust(d_dietoverlap, method = "complete")

dend <- as.dendrogram(guilds)

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

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

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


# Filter NEFSC food habits data with predator list
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))


##############################################################################
# Get prey list from NEFSC and NEAMAP

preycount  <- fh.nefsc.pisc.pisccomplete %>%
  #group_by(year, season, pdcomnam, pynam) %>%
  group_by(pdcomnam, pynam) %>%
  summarise(count = n()) %>%
  #arrange(desc(count))
  pivot_wider(names_from = pdcomnam, values_from = count)


gencomlist <- allfh %>%
  select(pynam, pycomnam2, gencom2) %>%
  distinct()

NEFSCblueprey <- preycount %>%
  #filter(BLUEFISH > 9) %>%
  filter(!pynam %in% c("EMPTY", "BLOWN",
                       "FISH", "OSTEICHTHYES",
                       "ANIMAL REMAINS",
                       "FISH SCALES")) %>%
  #filter(!str_detect(pynam, "SHRIMP|CRAB")) %>%
  left_join(gencomlist) %>%
  filter(!gencom2 %in% c("ARTHROPODA", "ANNELIDA",
                         "CNIDARIA", "UROCHORDATA",
                         "ECHINODERMATA", "WORMS",
                         "BRACHIOPODA", "COMB JELLIES",
                         "BRYOZOA", "SPONGES",
                         "MISCELLANEOUS", "OTHER")) %>%
  arrange(desc(BLUEFISH))

# March 2023, formally add NEAMAP to prey decisions
NEAMAPblueprey <- read.csv(here("fhdat/Full Prey List_Common Names.csv")) %>%
  #filter(BLUEFISH > 9) %>%
  filter(!SCIENTIFIC.NAME %in% c("Actinopterygii", "fish scales",
                                 "Decapoda (megalope)", 
                                 "unidentified material",
                                 "Plantae",
                                 "unidentified animal"))

NEAMAPprey <- NEAMAPblueprey %>%
  dplyr::select(COMMON.NAME, SCIENTIFIC.NAME, BLUEFISH) %>%
  dplyr::filter(!is.na(BLUEFISH)) %>%
  dplyr::mutate(pynam2 = tolower(SCIENTIFIC.NAME),
                pynam2 = stringr::str_replace(pynam2, "spp.", "sp")) %>%
  dplyr::rename(NEAMAP = BLUEFISH)


NEFSCprey <- NEFSCblueprey %>%
  dplyr::select(pycomnam2, pynam, BLUEFISH) %>%
  dplyr::filter(!is.na(BLUEFISH)) %>%
  dplyr::mutate(pynam2 = tolower(pynam)) %>%
  dplyr::rename(NEFSC = BLUEFISH)

# new criteria March 2023, >20 observations NEAMAP+NEFSC, but keep mackerel
# removes the flatfish order (too broad) and unid Urophycis previously in NEAMAP
blueprey <- NEFSCprey %>% 
  dplyr::full_join(NEAMAPprey) %>%
  dplyr::mutate(NEAMAP = ifelse(is.na(NEAMAP), 0, NEAMAP),
                NEFSC = ifelse(is.na(NEFSC), 0, NEFSC),
                total = NEFSC + NEAMAP,
                PREY = ifelse(is.na(SCIENTIFIC.NAME), pynam, SCIENTIFIC.NAME),
                COMMON = ifelse(is.na(COMMON.NAME), pycomnam2, COMMON.NAME),
                pynam = ifelse(is.na(pynam), toupper(pynam2), pynam)) %>%
  dplyr::arrange(desc(total)) %>%
  dplyr::filter(total>20 | pynam=="SCOMBER SCOMBRUS") %>% # >20 leaves out mackerel
  dplyr::mutate(COMMON = case_when(pynam=="ILLEX SP" ~ "Shortfin squids",
                                   pynam2=="teuthida" ~ "Unidentified squids",
                                   TRUE ~ COMMON)) %>%
  dplyr::mutate(PREY = stringr::str_to_sentence(PREY),
                COMMON = stringr::str_to_sentence(COMMON))


fh.nefsc.pisc.pisccomplete.blueprey <- fh.nefsc.pisc.pisccomplete %>%
  mutate(blueprey = case_when(pynam %in% blueprey$pynam ~ "blueprey",
                              TRUE ~ "othprey"))

###############################################################################
# Make the NEFSC dataset aggregating prey based on prey list

bluepyall_stn <- fh.nefsc.pisc.pisccomplete.blueprey %>%
  #create id linking cruise6_station
  #create season_ng spring and fall Spring=Jan-May, Fall=June-Dec
  mutate(id = paste0(cruise6, "_", station),
         year = as.numeric(year),
         month = as.numeric(month),
         season_ng = case_when(month <= 6 ~ "SPRING",
                               month >= 7 ~ "FALL",
                               TRUE ~ as.character(NA))
  ) %>%
  dplyr::select(year, season_ng, id, stratum,
                pynam, pyamtw, pywgti, pyvoli, blueprey, 
                pdcomnam, pdid, pdlen, pdsvol, pdswgt, 
                beglat, beglon, declat, declon, 
                bottemp, surftemp, setdepth) %>%
  group_by(id) %>%
  #mean blueprey g per stomach per tow: sum all blueprey g/n stomachs in tow
  mutate(bluepywt = case_when(blueprey == "blueprey" ~ pyamtw,
                              TRUE ~ 0.0),
         bluepynam = case_when(blueprey == "blueprey" ~ pynam,
                               TRUE ~ NA_character_)) 

# Optional: save at prey disaggregated stage for paper
#saveRDS(bluepyall_stn, here("fhdat/bluepyall_stn.rds"))

# Now get station data in one line
stndat <- bluepyall_stn %>%
  dplyr::select(year, season_ng, id, 
                beglat, beglon, declat, declon, 
                bottemp, surftemp, setdepth) %>%
  distinct()

#pisc stomachs in tow count pdid for each pred and sum
piscstom <- bluepyall_stn %>%
  group_by(id, pdcomnam) %>%
  summarise(nstompd = n_distinct(pdid)) %>%
  group_by(id) %>%
  summarise(nstomtot = sum(nstompd))

#mean and var pred length per tow
pisclen <- bluepyall_stn %>%
  summarise(meanpisclen = mean(pdlen),
            varpisclen = var(pdlen))

# Aggregated prey at station level with predator covariates
bluepyagg_stn <- bluepyall_stn %>%
  summarise(sumbluepywt = sum(bluepywt),
            nbluepysp = n_distinct(bluepynam, na.rm = T),
            npreysp = n_distinct(pynam),
            npiscsp = n_distinct(pdcomnam)) %>%
  left_join(piscstom) %>%
  mutate(meanbluepywt = sumbluepywt/nstomtot) %>%
  left_join(pisclen) %>%
  left_join(stndat)

# save at same stage as before, writing over old file
#saveRDS(bluepyagg_stn, here("fhdat/bluepyagg_stn.rds"))

# current dataset, fix declon, add vessel, rename NEFSC
#nefsc_bluepyagg_stn <- readRDS(here("fhdat/bluepyagg_stn.rds")) %>%
nefsc_bluepyagg_stn <- bluepyagg_stn %>%
  mutate(declon = -declon,
         vessel = case_when(year<2009 ~ "AL",
                            year>=2009 ~ "HB", 
                            TRUE ~ as.character(NA)))

##############################################################################
# Add NEAMAP to make full aggregated stomach dataset

# Read in NEAMAP updated input from Jim Gartland, reformat with same names
neamap_bluepreyagg_stn <- read_csv(here("fhdat/NEAMAP_Mean stomach weights_Bluefish Prey_Oct2023.csv")) %>%  
  mutate(vessel = "NEAMAP") %>%
  rename(id = station,
         sumbluepywt = sumbluepreywt,
         nbluepysp = nbfpreyspp,
         #npreysp = ,
         npiscsp = npiscspp,
         nstomtot = nstomtot, 
         meanbluepywt = meanbluepreywt,
         meanpisclen = meanpisclen.simple, 
         #varpisclen = ,
         season_ng = season,
         declat  = lat,
         declon = lon,
         bottemp = bWT,
         #surftemp = , 
         setdepth = depthm) 

# Read in NEAMAP 2023 input from Jim Gartland, reformat with same names
neamap_bluepreyagg_stn23 <- read_csv(here("fhdat/NEAMAP_Mean stomach weights_Bluefish Prey_Oct2024.csv")) %>%  
  mutate(vessel = "NEAMAP") %>%
  rename(id = station,
         sumbluepywt = sumbluepreywt,
         nbluepysp = nbfpreyspp,
         #npreysp = ,
         npiscsp = npiscspp,
         nstomtot = nstomtot, 
         meanbluepywt = meanbluepreywt,
         meanpisclen = meanpisclen.simple, 
         #varpisclen = ,
         season_ng = season,
         declat  = lat,
         declon = lon,
         bottemp = bWT,
         #surftemp = , 
         setdepth = depthm) 

neamap_bluepreyagg_stn <- dplyr::bind_rows(neamap_bluepreyagg_stn, neamap_bluepreyagg_stn23)

# combine NEAMAP and NEFSC
bluepyagg_stn_all <-  nefsc_bluepyagg_stn %>%
  bind_rows(neamap_bluepreyagg_stn) 

# Save before SST integration step
#saveRDS(bluepyagg_stn_all, here("fhdat/bluepyagg_stn_all.rds"))

###############################################################################
# Add SST into NEAMAP and reintegrate into full dataset

# Read back in if needed for SST
#bluepyagg_stn_all <- readRDS(here("fhdat/bluepyagg_stn_all.rds"))

NEFSCstations <- allfh %>%
  dplyr::mutate(id = paste0(cruise6, "_", station),
                year = as.numeric(year),
                month = as.numeric(month),
                day = as.numeric(day),
                declon = -declon) %>%
  dplyr::select(id, year, month, day, declat, declon) %>%
  dplyr::distinct()

# Need NEAMAP SST update! This is the old file
NEAMAPstationSST <- read.csv(here("fhdat/NEAMAP SST_2007_2022.csv"))

NEAMAPstationSST23 <- read.csv(here("fhdat/NEAMAP SST_2023.csv"))

NEAMAPstationSST <- dplyr::bind_rows(NEAMAPstationSST, NEAMAPstationSST23)

NEAMAPstations <- NEAMAPstationSST %>%
  dplyr::mutate(id = station,
                year = as.numeric(year),
                month = as.numeric(month),
                day = as.numeric(day)) %>%
  dplyr::select(id, year, month, day) %>%
  dplyr::distinct()

# remake diethauls
diethauls <- bluepyagg_stn_all %>%
  dplyr::select(id, declat, declon)

NEFSCstations <- dplyr::select(NEFSCstations, c(-declat, -declon))

Allstations <- bind_rows(NEFSCstations, NEAMAPstations)

#station id, lat lon, year month day
diethauls <- left_join(diethauls, Allstations)

#add year month day to diet data
bluepyagg_stn_all <- left_join(bluepyagg_stn_all, diethauls)

# add NEAMAP SST to surftemp field
NEAMAPidSST <- NEAMAPstationSST %>%
  mutate(id = station) %>%
  dplyr::select(id, SST)

bluepyagg_stn_all <- left_join(bluepyagg_stn_all, NEAMAPidSST, by="id") %>%
  mutate(surftemp = coalesce(surftemp, SST)) %>%
  dplyr::select(-SST)

# save merged dataset with day, month, and NEAMAP surftemp, same name
#saveRDS(bluepyagg_stn_all, here("fhdat/bluepyagg_stn_all.rds"))

###############################################################################
#Now match stations to OISST

#make an SST dataframe for 2023! Add to saved sst_data in data-raw/gridded

library(sf)
library(raster)
library(terra)
library(nngeo)

# Bastille function from https://github.com/kimberly-bastille/ecopull/blob/main/R/utils.R

nc_to_raster <- function(nc,
                         varname,
                         extent = c(0, 360, -90, 90),
                         crop = raster::extent(280, 300, 30, 50),
                         show_images = FALSE) {
  
  message("Reading .nc as brick...")
  
  r <- raster::brick(nc, varname = varname)
  
  message("Setting CRS...")
  raster::crs(r) <- "+proj=longlat +lat_1=35 +lat_2=45 +lat_0=40 +lon_0=-77 +x_0=0 +y_0=0 +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"
  
  # not sure if this is necessary?
  raster::extent(r) <- raster::extent(extent)
  
  if(show_images){
    par(mfrow = c(1,2))
    raster::plot(r, 1, sub = "Full dataset")
  }
  
  message("Cropping data...")
  ne_data <- raster::crop(r, crop)
  #ne_data <- raster::rotate(ne_data) add here for future pulls
  
  if(show_images){
    raster::plot(ne_data, 1, sub = "Cropped dataset")
    par(mfrow = c(1,1))
  }
  
  message("Done!")
  
  return(ne_data)
}

# function to convert to dataframe based on
# https://towardsdatascience.com/transforming-spatial-data-to-tabular-data-in-r-4dab139f311f

raster_to_sstdf <- function(brick,
                            rotate=TRUE){
  
  if(rotate) brick_r <- raster::rotate(brick)
  brick_r <- raster::crop(brick_r, raster::extent(-77,-65,35,45))
  sstdf <- as.data.frame(raster::rasterToPoints(brick_r, spatial = TRUE))
  sstdf <- sstdf %>%
    dplyr::rename(Lon = x,
                  Lat = y) %>%
    tidyr::pivot_longer(cols = starts_with("X"),
                        names_to = c("year", "month", "day"),
                        names_prefix = "X",
                        names_sep = "\\.",
                        values_to = "sst",
    )
  return(sstdf)
}

# pull the OISST data as raster brick, modified from 
# https://github.com/kimberly-bastille/ecopull/blob/main/.github/workflows/pull_satellite_data.yml

varname <- "sst"

# 1985-2022 previously pulled, processed and stored. add 2023.
# add 1981-1984 to extend back in time. No OISST before 1981.
# 1981 is only Sept-Dec so don't use

years <- 2023 #1982:1984 # 2022
for(i in years) {
  name <- paste0(i, ".nc")
  dir.create(here::here("data-raw","gridded", "sst_data"), recursive = TRUE)
  filename <- here::here("data-raw","gridded", "sst_data", paste0("test_", i, ".grd"))
  url <- paste0("https://downloads.psl.noaa.gov/Datasets/noaa.oisst.v2.highres/sst.day.mean.", i, ".nc")
  options(timeout = max(300, getOption("timeout")))
  download.file(url, destfile = name)
  
  text <- knitr::knit_expand(text = "test_{{year}} <- nc_to_raster(nc = name, varname = varname)
                                     raster::writeRaster(test_{{year}}, filename = filename, overwrite=TRUE)",
                             year = i)
  print(text)
  try(eval(parse(text = text)))
  unlink(name) # remove nc file to save space
  print(paste("finished",i))
}


# convert raster to dataframe
#years <- 2022
for(i in years) {
  name <- get(paste0("test_",i))
  filename <- here::here("data-raw","gridded", "sst_data", paste0("sst", i, ".rds"))
  text <- knitr::knit_expand(text = "sst{{year}} <- raster_to_sstdf(brick = name)
                                     saveRDS(sst{{year}}, filename)",
                             year = i)
  print(text)
  try(eval(parse(text = text)))
}
#read in diet data with month day fields

#bluepyagg_stn_all <- readRDS(here("fhdat/bluepyagg_stn_all.rds"))

stations <- bluepyagg_stn_all %>%
  dplyr::mutate(day = str_pad(day, 2, pad='0'),
                month = str_pad(month, 2, pad='0'),
                yrmody = as.numeric(paste0(year, month, day))) %>%
  dplyr::select(id, declon, declat, year, yrmody) %>%
  na.omit() %>%
  sf::st_as_sf(coords=c("declon","declat"), crs=4326, remove=FALSE) 



#list of SST dataframes
SSTdfs <- list.files(here("data-raw/gridded/sst_data/"), pattern = "*.rds")

dietstn_OISST <- tibble()


for(df in SSTdfs){
  sstdf <- readRDS(paste0(here("data-raw/gridded/sst_data/", df)))
  
  # keep only bluefish dates in SST year
  stationsyr <- stations %>%
    filter(year == unique(sstdf$year))
  
  # keep only sst days in bluefish dataset
  sstdf_survdays <- sstdf %>%
    dplyr::mutate(yrmody = as.numeric(paste0(year, month, day)) )%>%
    dplyr::filter(yrmody %in% unique(stationsyr$yrmody)) %>%
    dplyr::mutate(year = as.numeric(year),
                  month = as.numeric(month),
                  day = as.numeric(day),
                  declon = Lon,
                  declat = Lat) %>%
    dplyr::select(-Lon, -Lat) %>%
    sf::st_as_sf(coords=c("declon","declat"), crs=4326, remove=FALSE)  
  
  # now join by nearest neighbor and date
  
  #https://stackoverflow.com/questions/71959927/spatial-join-two-data-frames-by-nearest-feature-and-date-in-r      
  
  yrdietOISST <- do.call('rbind', lapply(split(stationsyr, 1:nrow(stationsyr)), function(x) {
    sf::st_join(x, sstdf_survdays[sstdf_survdays$yrmody == unique(x$yrmody),],
                #join = st_nearest_feature
                join = st_nn, k = 1, progress = FALSE
    )
  }))
  
  #   #datatable solution--works but doesnt seem faster?
  #    df1 <- data.table(stationsyr)
  #   
  #  .nearest_samedate <- function(x) {
  #    st_join(st_as_sf(x), sstdf_survdays[sstdf_survdays$yrmody == unique(x$yrmody),], join = st_nearest_feature)
  #  }
  # # 
  #  yrdietOISST <- df1[, .nearest_samedate(.SD), by = list(1:nrow(df1))]
  
  dietstn_OISST <- rbind(dietstn_OISST, yrdietOISST)
  
}

# takes >10 minutes to run, save and tack on next years rather than full merge?
saveRDS(dietstn_OISST, here("data-raw/dietstn_OISST_1982_2023.rds"))

# Now join with OISST dataset

#bluepyagg_stn_all <- readRDS(here("fhdat/bluepyagg_stn_all.rds"))
#dietstn_OISST <- readRDS(here("data-raw/dietstn_OISST.rds"))


dietstn_OISST_merge <- dietstn_OISST %>%
  dplyr::rename(declon = declon.x,
                declat = declat.x,
                year = year.x,
                oisst = sst) %>%
  dplyr::select(id, oisst) %>%
  sf::st_drop_geometry()

bluepyagg_stn_all_OISST <- left_join(bluepyagg_stn_all, dietstn_OISST_merge)

#saveRDS(bluepyagg_stn_all_OISST, here("fhdat/bluepyagg_stn_all_OISST_1982-2022.rds"))
saveRDS(bluepyagg_stn_all_OISST, here("fhdat/bluepyagg_stn_all_OISST_1982-2023.rds"))

For the 2024 report, VAST models were run using the script VASTunivariate_seasonalforageindex_operational.R in the folder VASTscripts: https://github.com/NOAA-EDAB/forageindex/blob/main/VASTscripts/VASTunivariate_seasonalforageindex_operational.R

That previous script fails for me under VAST v. ‘3.11.2’ and FishStatsUtils v. ‘2.13.1’. I also had to upgrade R to 4.4.1, all of this to run WHAM.

So, for the 2025 report I have used this script to get our indices in on time: https://github.com/NOAA-EDAB/forageindex/blob/main/VASTscripts/VASTunivariate_seasonalforageindex_SOE_EPUsonly.R

# This is the VAST code used in Gaichas et al 2023, but with data years 
# extended from 1973-2023
# AND NO CUSTOM EXTRAPOLATION GRID SO NO 3-nmile INSHORE STRATA
# VAST attempt 2 univariate model as a script
# modified from https://github.com/James-Thorson-NOAA/VAST/wiki/Index-standardization

# Load packages
library(here)
library(dplyr)
library(VAST)

#Read in data, separate spring and fall, and rename columns for VAST:

# this dataset created in SSTmethods.Rmd

bluepyagg_stn <- readRDS(here::here("fhdat/bluepyagg_stn_all_OISST_1982-2023.rds"))

# 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, #leaves out everything before 1982
         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, #leaves out everything before 1982
         sstfill
         ) %>%
  na.omit() %>%
  as.data.frame()


# Make settings (turning off bias.correct to save time for example)
# NEFSC strata limits https://github.com/James-Thorson-NOAA/VAST/issues/302

# use only MAB, GB, GOM, SS EPUs 
# leave out south of Cape Hatteras at Elizabeth's suggestion
# could also leave out SS?
# CHECK if these EPUs match what we use in SOE

bfinshore <- c(3020, 3050, 3080, 3110, 3140, 3170, 3200, 3230, 
              3260, 3290, 3320, 3350, 3380, 3410, 3440, 3450, 3460)

bfoffshore <- c(1010, 1730, 1690, 1650, 1050, 1060, 1090, 1100, 1250, 1200, 1190, 1610)

MAB <- c(1010:1080, 1100:1120, 1600:1750, 3010:3450, 3470, 3500, 3510)
GB  <- c(1090, 1130:1210, 1230, 1250, 3460, 3480, 3490, 3520:3550)
GOM <- c(1220, 1240, 1260:1290, 1360:1400, 3560:3830)
SS  <- c(1300:1352, 3840:3990)

# using VAST built in grid because custom one is breaking Oct 2024
MAB2 <- FishStatsUtils::northwest_atlantic_grid %>% 
  dplyr::filter(stratum_number %in% MAB) %>%
  dplyr::select(stratum_number) %>%
  dplyr::distinct()

# Georges Bank EPU
GB2 <- FishStatsUtils::northwest_atlantic_grid %>% 
  dplyr::filter(stratum_number %in% GB) %>%
  dplyr::select(stratum_number) %>%
  dplyr::distinct()

# gulf of maine EPU (for SOE)
GOM2 <- FishStatsUtils::northwest_atlantic_grid %>% 
  dplyr::filter(stratum_number %in% GOM) %>%
  dplyr::select(stratum_number) %>%
  dplyr::distinct()

# scotian shelf EPU (for SOE)
SS2 <- FishStatsUtils::northwest_atlantic_grid %>%  
  dplyr::filter(stratum_number %in% SS) %>%
  dplyr::select(stratum_number) %>%
  dplyr::distinct()

# needed to cover the whole northwest atlantic grid--lets try without
 allother2 <- FishStatsUtils::northwest_atlantic_grid %>%
   dplyr::filter(!stratum_number %in% c(MAB, GB, GOM, SS)) %>%
   dplyr::select(stratum_number) %>%
   dplyr::distinct()

# all epus
allEPU2 <- FishStatsUtils::northwest_atlantic_grid %>%  
  dplyr::filter(stratum_number %in% c(MAB, GB, GOM, SS)) %>%
  dplyr::select(stratum_number) %>%
  dplyr::distinct()

# # configs same FieldConfig as below but formatted differently
# FieldConfig <- c(
#   "Omega1"   = 0,   # number of spatial variation factors (0, 1, AR1)
#   "Epsilon1" = 0,   # number of spatio-temporal factors
#   "Omega2"   = 0, 
#   "Epsilon2" = 0
# ) 

# default configs, not really specified anyway
FieldConfig = matrix( "IID", ncol=2, nrow=3, 
                      dimnames=list(c("Omega","Epsilon","Beta"),c("Component_1","Component_2")))


RhoConfig <- c(
  "Beta1" = 0,      # temporal structure on years (intercepts) 
  "Beta2" = 0, 
  "Epsilon1" = 0,   # temporal structure on spatio-temporal variation
  "Epsilon2" = 0
) 
# 0 off (fixed effects)
# 1 independent
# 2 random walk
# 3 constant among years (fixed effect)
# 4 AR1

OverdispersionConfig    <- c("eta1"=0, "eta2"=0)
# eta1 = vessel effects on prey encounter rate
# eta2 = vessel effects on prey weight

strata.limits <- as.list(c("AllEPU" = allEPU2, 
                           "MAB" = MAB2,
                           "GB" = GB2,
                           "GOM" = GOM2,
                           "SS" = SS2,
                           "allother" = allother2
                           ))

settings = make_settings( n_x = 500, 
                          Region = "northwest_atlantic",
                          #Version = "VAST_v14_0_1", #needed to prevent error from newer dev version number
                          #strata.limits = list('All_areas' = 1:1e5), full area
                          strata.limits = strata.limits,
                          purpose = "index2", 
                          bias.correct = TRUE,
                          #use_anisotropy = FALSE,
                          #fine_scale = FALSE,
                          #FieldConfig = FieldConfig,
                          #RhoConfig = RhoConfig,
                          OverdispersionConfig = OverdispersionConfig
                          )


#New_Extrapolation_List <- readRDS(here::here("spatialdat/CustomExtrapolationList.rds"))

# select dataset and set directory for output

#########################################################
# Run model fall

season <- c("fall_500_lennosst_EPUonly_biascorrect")

working_dir <- here::here(sprintf("SOEpyindex/allagg_%s/", season))

if(!dir.exists(working_dir)) {
  dir.create(working_dir)
}

# subset for faster testing
#bluepyagg_stn_fall <- bluepyagg_stn_fall %>% filter(Year<1990)                       

fit <- fit_model(
  settings = settings, 
  #extrapolation_list = New_Extrapolation_List,
  Lat_i = bluepyagg_stn_fall$Lat, 
  Lon_i = bluepyagg_stn_fall$Lon, 
  t_i = bluepyagg_stn_fall$Year, 
  b_i = as_units(bluepyagg_stn_fall[,'Catch_g'], 'g'),
  a_i = rep(1, nrow(bluepyagg_stn_fall)),
  v_i = bluepyagg_stn_fall$Vessel,
  Q_ik = as.matrix(bluepyagg_stn_fall[,c("npiscsp", 
                                         "meanpisclen", 
                                         "sstfill"
                                         )]),
  #Use_REML = TRUE,
  working_dir = paste0(working_dir, "/"))

saveRDS(fit, file = paste0(working_dir, "/fit.rds"))

# Plot results
plot( fit,
      working_dir = paste0(working_dir, "/"))

######################################################
# Run model spring

season <- c("spring_500_lennosst_EPUonly_biascorrect")

working_dir <- here::here(sprintf("SOEpyindex/allagg_%s/", season))

if(!dir.exists(working_dir)) {
  dir.create(working_dir)
}                         
  
# subset for faster testing
#bluepyagg_stn_spring <- bluepyagg_stn_spring %>% filter(Year<1990) 

fit <- fit_model( settings = settings,  
                 #extrapolation_list = New_Extrapolation_List,
                 Lat_i = bluepyagg_stn_spring[,'Lat'], 
                 Lon_i = bluepyagg_stn_spring[,'Lon'], 
                 t_i = bluepyagg_stn_spring[,'Year'], 
                 b_i = as_units(bluepyagg_stn_spring[,'Catch_g'], 'g'), 
                 a_i = rep(1, nrow(bluepyagg_stn_spring)),
                 v_i = bluepyagg_stn_spring$Vessel,
                 Q_ik = as.matrix(bluepyagg_stn_spring[,c("npiscsp", 
                                                          "meanpisclen", 
                                                          "sstfill"
                                                          )]),
                # Use_REML = TRUE,
                 working_dir = paste0(working_dir, "/"))

saveRDS(fit, file = paste0(working_dir, "/fit.rds"))

# Plot results
plot( fit,
      working_dir = paste0(working_dir, "/")) 

Model output was saved in the folder SOEpyindex for the time period submitted in 2024 plus the updated data: 1982-2023.

The script to create the SOE forage indices from the VAST output is in the SOEpyindex folder: https://github.com/NOAA-EDAB/forageindex/blob/main/SOEpyindex/SOE-VASTForageIndices.R

# create csv for ecodata input
# aim for similar structure to other ecodata datasets


library(dplyr)
library(ggplot2)
library(tidyr)

SOEinputs <- function(infile, season, stratlook, outfile) {
  
  splitoutput <- read.csv(infile)
  
  stratlook <- stratlook
  
  forageindex <- splitoutput %>%
    left_join(stratlook) %>%
    dplyr::select(Time, 
                  EPU = Region, 
                  "Forage Fish Biomass Estimate" = Estimate, 
                  "Forage Fish Biomass Estimate SE" = Std..Error.for.Estimate) %>%
    tidyr::pivot_longer(c("Forage Fish Biomass Estimate", "Forage Fish Biomass Estimate SE"), 
                        names_to = "Var", values_to = "Value") %>%
    dplyr::filter(EPU %in% c("MAB", "GB", "GOM", "AllEPU")) %>%
    dplyr::mutate(Units = "relative grams per stomach") %>%
    dplyr::select(Time, Var, Value, EPU, Units)
  
  forageindex$Var <- stringr::str_c(season, forageindex$Var, sep = " ")
  
  #readr::write_csv(forageindex, outfile)
  saveRDS(forageindex, outfile)
  
} 

# warning, hardcoded. obviously
stratlook_ALLsplit <- 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"))

stratlook_EPUonly <- data.frame(Stratum = c("Stratum_1",
                                            "Stratum_2",
                                            "Stratum_3",
                                            "Stratum_4",
                                            "Stratum_5",
                                            "Stratum_6"),
                                Region = c("AllEPU", 
                                           "MAB",
                                           "GB",
                                           "GOM",
                                           "SS",
                                           "allother"))



# make data files
SOEinputs(infile = "SOEpyindex/1982-2022/allagg_fall_500_lennosst_ALLsplit_biascorrect/Index.csv",
          season = "Fall", 
          stratlook = stratlook_EPUonly, 
           outfile = "SOEpyindex/fallforageindex.rds")

SOEinputs(infile = "SOEpyindex/1982-2022/allagg_spring_500_lennosst_ALLsplit_biascorrect/Index.csv",
          season = "Spring", 
          stratlook = stratlook_EPUonly,
           outfile = "SOEpyindex/springforageindex.rds")

# SOEinputs(infile = "pyindex/allagg_annual_500_lennosst_ALLsplit_biascorrect/Index.csv",
#           season = "Annual",
#           stratlook = stratlook_ALLsplit,
#            outfile = "toSOE/annualforageindex.rds")



# test plot
# foragewide <- forageindex %>%
#   pivot_wider(names_from = Var, values_from = Value)
# 
# 
# ggplot(foragewide, aes(x=Time, y=`Forage Fish Biomass Estimate`, colour = EPU)) +
#   geom_errorbar(aes(ymin=`Forage Fish Biomass Estimate`+`Forage Fish Biomass Estimate SE`, 
#                     ymax=`Forage Fish Biomass Estimate`-`Forage Fish Biomass Estimate SE`))+
#   geom_point()+
#   geom_line()

Code blocks below setup the index comparisons in the remainder of this document, where I compare only the 1982-2022 (2024 SOE) and 1982-2023 (2025 SOE) indices.

# compare full custom strata set to to EPU only strata

stratlook <- data.frame(Stratum = c("ALLsplit_Stratum_1",
                                    "ALLsplit_Stratum_2",
                                    "ALLsplit_Stratum_3",
                                    "ALLsplit_Stratum_4",
                                    "ALLsplit_Stratum_5",
                                    "ALLsplit_Stratum_6",
                                    "ALLsplit_Stratum_7",
                                    "ALLsplit_Stratum_8",
                                    "ALLsplit_Stratum_9",
                                    "ALLsplit_Stratum_10",
                                    "ALLsplit_Stratum_11",
                                    "ALLsplit_Stratum_12",
                                    "ALLsplit_Stratum_13",
                                    "ALLsplit_Stratum_14",
                                    "ALLsplit_Stratum_15"),
                        Region  = c("AllEPU", 
                                    "MABGB", 
                                    "MABGBstate", 
                                    "MABGBfed", 
                                    "MAB",
                                    "GB",
                                    "GOM",
                                    "bfall",
                                    "bfin",
                                    "bfoff",
                                    "MABGBalbinshore",
                                    "MABGBothoffshore",
                                    "albbfin",
                                    "albbfall",
                                    "allother"))
# new run has fewer strata
stratlook2 <- data.frame(Stratum = c("EPUonly_Stratum_1",
                                     "EPUonly_Stratum_2",
                                     "EPUonly_Stratum_3",
                                     "EPUonly_Stratum_4",
                                     "EPUonly_Stratum_5",
                                     "EPUonly_Stratum_6"),
                         Region = c("AllEPU", 
                                    "MAB",
                                    "GB",
                                    "GOM",
                                    "SS",
                                    "allother"))

stratlookall <- dplyr::bind_rows(stratlook, stratlook2) 

# from each output folder in pyindex, 
outdir <- here("SOEpyindex")
moddirs <- list.dirs(outdir)
# only get pred sensitivities and original full models
fullfall <- grep("/allagg_fall_", moddirs, value = TRUE)
fullspring <- grep("/allagg_spring_", moddirs, value = TRUE)
moddirs <- c(fullfall, fullspring)
# get original published version 1985-2021
#outdir1 <- here("pyindex")
#moddirs1 <- list.dirs(outdir1)
#RTfall <- grep("pyindex/2022/allagg_fall_", moddirs1, value = TRUE)
#RTspring <- grep("pyindex/2022/allagg_spring_", moddirs1, value = TRUE)
#paperfall <- grep("pyindex/allagg_fall_", moddirs1, value = TRUE)
#paperspring <- grep("pyindex/allagg_spring_", moddirs1, value = TRUE)
#moddirs1 <- c(RTfall, RTspring, paperfall, paperspring)
# keep folder name
modnames <- list.dirs(outdir, full.names = FALSE)
#modnames1 <- list.dirs(outdir1, 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)]
  modnamepre <- modpath[length(modpath)-1]
  
  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 = paste0(modnamepre,"_", 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)
# }


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

moddirsall <- c(moddirs)#, moddirs1)

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

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

3 Check for convergence, etc

modtable <- function(moddirs){
  
  # apply getmodinfo function to inout directories
  modcompare <- purrr::map_dfr(moddirs, getmodinfo)
  
  modselect <- modcompare |>
    dplyr::mutate(season = dplyr::case_when(stringr::str_detect(modname, "_fall_") ~ "Fall",
                                            stringr::str_detect(modname, "spring") ~ "Spring",
                                            stringr::str_detect(modname, "_all_") ~ "Annual",
                                            TRUE ~ as.character(NA))) |>
    dplyr::mutate(converged2 = dplyr::case_when(stringr::str_detect(converged, "no evidence") ~ "likely",
                                                stringr::str_detect(converged, "is likely not") ~ "unlikely",
                                                TRUE ~ as.character(NA))) |>
    dplyr::mutate(timeseries = stringr::str_extract(modname, "[^_]+")) |>
    #dplyr::mutate(modname = str_extract(modname, '(?<=allagg_).*')) |>
    dplyr::group_by(timeseries, season) |>
    #dplyr::mutate(deltaAIC = AIC-min(AIC)) |>
    dplyr::select(timeseries, modname, season, #deltaAIC, 
                  fixedcoeff,
                  randomcoeff, use_anisotropy, 
                  omega1, omega2, epsilon1, epsilon2, 
                  beta1, beta2, AIC, converged2) |>
    dplyr::arrange(timeseries, season, AIC)
  
  return(modselect)
}

modselect <- modtable(moddirs)

flextable::flextable(modselect %>%
                       dplyr::mutate(modname = str_extract(modname, '(?<=allagg_).*')) |>
                       dplyr::filter(timeseries %in% c("1982-2022", "1982-2023")) |>
                       dplyr::select(-c(use_anisotropy, 
         omega1, omega2, epsilon1, epsilon2,
         beta1, beta2))) %>%
  flextable::set_header_labels(timeseries = "Data",
                               modname = "Model name",
                               season = "Season",
                               #deltaAIC = "dAIC",
                               fixedcoeff = "N fixed",
                               randomcoeff = "N random",
                               converged2 = "Convergence") |>
  flextable::set_caption("Convergence check for 2025 SOE submission relative to 2024.") %>%
  flextable::fontsize(size = 9, part = "all") %>%
  flextable::colformat_double(digits = 2) |>
  flextable::set_table_properties(layout = "autofit", width = 1)

``` ```{=html}

(\#tab:unnamed-chunk-4)Convergence check for 2025 SOE submission relative to 2024.

Data	Model name	Season	N fixed	N random	AIC	Convergence
1982-2022	fall_500_lennosst_ALLsplit_biascorrect	Fall	97	46,368	52,682.11	likely
1982-2022	spring_500_lennosst_ALLsplit_biascorrect	Spring	97	46,368	53,302.56	likely
1982-2023	fall_500_lennosst_EPUonly_biascorrect	Fall	99	47,472	53,998.53	likely
1982-2023	spring_500_lennosst_EPUonly_biascorrect	Spring	99	47,558	54,132.73	likely

# # lets only look at indices for converged models
# modcompare_conv <- modselect |>
#   dplyr::ungroup() |>
#   dplyr::filter(converged2 == "likely") |>
#   dplyr::select(modname) |>
#   as.vector() |>
#   unname() |>
#   unlist()
# 
# moddirs_conv <-  moddirs[grepl(sprintf("\\.*(%s)$", paste(modcompare_conv, collapse = '|')), moddirs)]

4 Trend comparisons

4.1 Methods change and adding one year

Lets compare last year’s 2024 SOE submission with the 2025 submission using the same model settings but without the custom extrapolation grid, and with one more year of data (2023).

This is fall trend.

splitoutput <- modcompareindex %>%
  dplyr::mutate(Season = modname |> map(str_split, pattern = "_") |> map_chr(c(1,3))) %>%
  dplyr::mutate(Data = modname |> map(str_split, pattern = "_") |> map_chr(c(1,1))) %>% 
  dplyr::mutate(modname = str_extract(modname, '(?<=allagg_).*')) |>
  tidyr::unite("Stratum", c(strata, Stratum) ) |>
  dplyr::left_join(stratlookall) %>%
  dplyr::filter(Region %in% c("GOM", "GB", "MAB", "AllEPU", "allother")) %>%
  dplyr::mutate(Type = ifelse(Region %in% c("GOM", "GB", "MAB"), "Ecoregion", "Coastwide"))

foragemax <- max(splitoutput$Estimate)


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


ggplot(foragetsmean |> filter(Season =="fall",
                             Data %in% c("1982-2022",
                                         "1982-2023")), 
       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)))

Figure 4.1: Trend comparison between fall forage indices using 1982-2022 (SOE 2024), and 1982-2023 (SOE 2025) with standard VAST grid.

This is spring trend.

ggplot(foragetsmean |> filter(Season =="spring",
                             Data %in% c("1982-2022",
                                         "1982-2023")), 
       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)))

Figure 4.2: Trend comparison between spring forage indices using 1982-2022 (SOE 2024), and 1982-2023 (SOE 2025) with standard VAST grid.

Trends look comparable. The difference looks like the time series mean changing with the addition of 2023. So far so good.

4.2 Stomach spatial coverage

Spatial coverage for stomachs in fall

newmods <- moddirs[stringr::str_detect(moddirs,"EPUonly")]

mapfall <- paste0(newmods[stringr::str_detect(newmods,"_fall_")], "/Data_by_year.png")

knitr::include_graphics(mapfall) 

Spatial coverage for stomachs in spring

mapspring <- paste0(newmods[stringr::str_detect(newmods,"_spring_")], "/Data_by_year.png")

knitr::include_graphics(mapspring) 

5 Scale comparisons

We don’t use scale much, but in case someone does here they are, scaled relative to the max observed across all datasets.

5.1 Methods change and adding one year

Same comparisons as above: fall

ggplot(splitoutput |> filter(Season =="fall",
                             Data %in% c("1982-2022",
                                         "1982-2023")), 
       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)))

Figure 5.1: Scale comparison between fall forage indices using 1982-2022 (SOE 2024), and 1982-2023 (SOE 2025) with standard VAST grid.

Spring

ggplot(splitoutput |> filter(Season =="spring",
                             Data %in% c("1982-2022",
                                         "1982-2023")), 
       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)))

Figure 5.2: Scale comparison between spring forage indices using 1982-2022 (SOE 2024), and 1982-2023 (SOE 2025) with standard VAST grid.

The spring comparison is what I was hoping for. I’m unclear on why the fall model is different in scale.

Plot them all on the reported scale

ggplot(splitoutput |> filter(#Season =="spring",
                             Data %in% c("1982-2022",
                                         "1982-2023")), 
       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~Season, scales = "free_y", ncol = 2,
              labeller = label_wrap_gen(multi_line=FALSE)) +
  #scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Forage biomass output by VAST")  +
  ggtitle("Output 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)))

6 Center of gravity exploration

Has forage shifted along the coast?

library(FishStatsUtils)

getcogVAST <- function(d.name){
  
  fit <- VAST::reload_model(readRDS(paste0(d.name,"/fit.rds")))
  
  dir.create(paste0(d.name,"/test"))
  
  cogout <- FishStatsUtils::plot_range_index(Sdreport = fit$parameter_estimates$SD, 
                                             Report = fit$Report, 
                                             TmbData = fit$data_list,
                                             year_labels = as.numeric(fit$year_labels),
                                             years_to_plot = fit$years_to_plot,
                                             Znames = colnames(fit$data_list$Z_xm),
                                             PlotDir = paste0(d.name,"/test")) #already have plots, will delete this directory
  
  saveRDS(cogout, paste0(d.name,"/cogout.rds"))
  
  unlink(paste0(d.name,"/test"), recursive=TRUE) #removes directory with unneeded plots
  
}

purrr::map(newmods, getcogVAST)

Fall center of gravity, forage index, significantly more North and East over time (same as last year)

cogoutfall <- readRDS(here::here("SOEpyindex/1982-2023/allagg_fall_500_lennosst_EPUonly_biascorrect/cogout.rds"))

cogdat <- as.data.frame(cogoutfall$COG_Table) |>
  dplyr::mutate(direction = ifelse(m==1, "Eastward", "Northward"))

ggplot2::ggplot(cogdat, ggplot2::aes(x = Year, y = COG_hat)) + 
  ggplot2::geom_point() + 
  ecodata::geom_gls() + 
  ggplot2::labs(y = "Center of gravity, km") +
  ggplot2::facet_wrap(~direction, scales = "free_y") + 
  ecodata::theme_facet()

Spring center of gravity, forage index, significantly northward with the addition of 2023 data (no significant trends last year)

cogoutspring <- readRDS(here::here("SOEpyindex/1982-2023/allagg_spring_500_lennosst_EPUonly_biascorrect/cogout.rds"))

cogdat <- as.data.frame(cogoutspring$COG_Table) |>
  dplyr::mutate(direction = ifelse(m==1, "Eastward", "Northward"))

ggplot2::ggplot(cogdat, ggplot2::aes(x = Year, y = COG_hat)) + 
  ggplot2::geom_point() + 
  ecodata::geom_gls() + 
  ggplot2::labs(y = "Center of gravity, km") +
  ggplot2::facet_wrap(~direction, scales = "free_y") + 
  ecodata::theme_facet()

update get function for ecodata to add cog to existing forage_index dataset

## Forage Index
raw.dir<- here::here("data-raw/")
#ann<- "annualforageindex - Sarah Gaichas - NOAA Federal.rds"
fal<- "fallforageindex - Sarah Gaichas - NOAA Federal.rds"
spr<- "springforageindex - Sarah Gaichas - NOAA Federal.rds"
falcog <- "fallforagecog.rds"
sprcog <- "springforagecog.rds"

get_forage_index <- function(save_clean = F){

  #annual<-readRDS(file.path(raw.dir, ann))
  fall<- readRDS(file.path(raw.dir, fal))
  spring<-readRDS(file.path(raw.dir, spr))
  fallcog <- readRDS(file.path(raw.dir, falcog))
  springcog <-readRDS(file.path(raw.dir, sprcog))

  forage_index<- rbind(fall, spring, fallcog, springcog)


  if (save_clean){
    usethis::use_data(forage_index, overwrite = T)
  } else {
    return(forage_index)
  }
}
get_forage_index(save_clean = T)

update plot function

plot_forage_index <- function(shadedRegion = NULL,
                              report="MidAtlantic",
                              varName = "index") {

  setup <- ecodata::plot_setup(shadedRegion = shadedRegion,
                               report=report)
  
  if (varName == "index") {
    
  if (report == "MidAtlantic") {
    filterEPUs <- c("MAB")
  } else {
    filterEPUs <- c("GB", "GOM")
  }

  fix<- ecodata::forage_index |>
    dplyr::filter(Var %in% c("Fall Forage Fish Biomass Estimate",
                             "Spring Forage Fish Biomass Estimate"),
                  EPU %in% filterEPUs) |>
    dplyr::group_by(EPU) |>
    dplyr::summarise(max = max(Value))

  p <- ecodata::forage_index |>
    dplyr::filter(Var %in% c("Fall Forage Fish Biomass Estimate",
                             "Fall Forage Fish Biomass Estimate SE",
                             "Spring Forage Fish Biomass Estimate",
                             "Spring Forage Fish Biomass Estimate SE"),
                  EPU %in% filterEPUs) |>
    dplyr::group_by(EPU) |>
    tidyr::separate(Var, into = c("Season", "A", "B", "C", "D", "Var")) |>
    dplyr::mutate(Var = tidyr::replace_na(Var, "Mean")) |> #,
                  #max = as.numeric(Value)) |>
    tidyr::pivot_wider(names_from = Var, values_from = Value) |>
    dplyr::left_join(fix) |>
    dplyr::mutate(#Value = Value/resca,
      Mean = as.numeric(Mean),
      #max = as.numeric(Value),
      Mean = Mean/max,
      SE = SE/max,
      Upper = Mean + SE,
      Lower = Mean - SE) |>
    ggplot2::ggplot(ggplot2::aes(x = Time, y = Mean, group = Season))+
    ggplot2::annotate("rect", fill = setup$shade.fill, alpha = setup$shade.alpha,
        xmin = setup$x.shade.min , xmax = setup$x.shade.max,
        ymin = -Inf, ymax = Inf) +
    ggplot2::geom_ribbon(ggplot2::aes(ymin = Lower, ymax = Upper, fill = Season), alpha = 0.5)+
    ggplot2::geom_point()+
    ggplot2::geom_line()+
    ggplot2::ggtitle("")+
    ggplot2::ylab(expression("Relative forage biomass"))+
    ggplot2::xlab(ggplot2::element_blank())+
    ggplot2::facet_wrap(.~EPU)+
    ecodata::geom_gls()+
    ecodata::theme_ts()+
    ecodata::theme_facet()+
    ecodata::theme_title()

    if (report == "NewEngland") {
      p <- p +
        ggplot2::theme(legend.position = "bottom",
                       legend.title = ggplot2::element_blank())

    }
  }
  
  if (varName == "cog"){
    
    p <- ecodata::forage_index |>
      dplyr::filter(Var %in% c("Fall Eastward Forage Fish Center of Gravity",
                               "Fall Eastward Forage Fish Center of Gravity SE",
                               "Fall Northward Forage Fish Center of Gravity",
                               "Fall Northward Forage Fish Center of Gravity SE",
                               "Spring Eastward Forage Fish Center of Gravity",
                               "Spring Eastward Forage Fish Center of Gravity SE",
                               "Spring Northward Forage Fish Center of Gravity",
                               "Spring Northward Forage Fish Center of Gravity SE")) |>
      tidyr::separate(Var, into = c("Season", "Direction", "A", "B", "C", "D", "E", "Var")) |>
      dplyr::mutate(Var = tidyr::replace_na(Var, "Mean")) |> 
      tidyr::pivot_wider(names_from = Var, values_from = Value) |>
      dplyr::mutate(Mean = as.numeric(Mean),
                    Upper = Mean + SE,
                    Lower = Mean - SE) |>
      ggplot2::ggplot(ggplot2::aes(x = Time, y = Mean, group = Season))+
      ggplot2::annotate("rect", fill = setup$shade.fill, alpha = setup$shade.alpha,
                        xmin = setup$x.shade.min , xmax = setup$x.shade.max,
                        ymin = -Inf, ymax = Inf) +
      ggplot2::geom_ribbon(ggplot2::aes(ymin = Lower, ymax = Upper, fill = Season), alpha = 0.3)+ #
      ggplot2::geom_point()+
      ggplot2::geom_line()+
      ggplot2::ggtitle("")+
      ggplot2::ylab(expression("Forage Center of Gravity, km"))+
      ggplot2::xlab(ggplot2::element_blank())+
      ggplot2::facet_wrap(~Direction, scales = "free_y")+ #Season
      ecodata::geom_gls()+
      ecodata::theme_ts()+
      ecodata::theme_facet()+
      ecodata::theme_title()

  }

    return(p)

}

attr(plot_forage_index,"report") <- c("MidAtlantic","NewEngland")
attr(plot_forage_index, "varName") <- c("index", "cog")

forage COG indices produced by SOE-VASTForageCOG.R

fallforagecog <- readRDS(here::here("SOEpyindex/fallforagecog.rds"))
springforagecog <- readRDS(here::here("SOEpyindex/springforagecog.rds"))

setup <- ecodata::plot_setup(shadedRegion = c(2014,2023), report = "MidAtlantic")

testplot <- dplyr::bind_rows(fallforagecog, springforagecog) |>
  
  dplyr::filter(Var %in% c("Fall Eastward Forage Fish Center of Gravity",
                               "Fall Eastward Forage Fish Center of Gravity SE",
                               "Fall Northward Forage Fish Center of Gravity",
                               "Fall Northward Forage Fish Center of Gravity SE",
                               "Spring Eastward Forage Fish Center of Gravity",
                               "Spring Eastward Forage Fish Center of Gravity SE",
                               "Spring Northward Forage Fish Center of Gravity",
                               "Spring Northward Forage Fish Center of Gravity SE")) |>
      tidyr::separate(Var, into = c("Season", "Direction", "A", "B", "C", "D", "E", "Var")) |>
      dplyr::mutate(Var = tidyr::replace_na(Var, "Mean")) |> 
      tidyr::pivot_wider(names_from = Var, values_from = Value) |>
      dplyr::mutate(Mean = as.numeric(Mean),
                    Upper = Mean + SE,
                    Lower = Mean - SE) |>
      ggplot2::ggplot(ggplot2::aes(x = Time, y = Mean, group = Season))+
      ggplot2::annotate("rect", fill = setup$shade.fill, alpha = setup$shade.alpha,
                        xmin = setup$x.shade.min , xmax = setup$x.shade.max,
                        ymin = -Inf, ymax = Inf) +
      ggplot2::geom_ribbon(ggplot2::aes(ymin = Lower, ymax = Upper, fill = Season), alpha = 0.3)+ #
      ggplot2::geom_point()+
      ggplot2::geom_line()+
      ggplot2::ggtitle("")+
      ggplot2::ylab(expression("Forage Center of Gravity, km"))+
      ggplot2::xlab(ggplot2::element_blank())+
      ggplot2::facet_wrap(~Direction, scales = "free_y")+ #Season
      ecodata::geom_gls()+
      ecodata::theme_ts()+
      ecodata::theme_facet()+
      ecodata::theme_title()

testplot