Having built an initial dataset of haul-specific mean bluefish prey weight per stomach for the updated piscivore list from both NEFSC and NEAMAP, and having completed initial model selection to see if spatial and spatio-temporal random effects make sense to estimate the prey index in VAST, now we’ll evaluate which catchability covariates are best supported by the data.

Following what Ng et al. (2021) did for herring, we apply a Poisson-link delta model to estimate expected prey mass per predator stomach. However, we use 500 knots, estimated by k-means clustering of the data, to define the spatial dimensions of each seasonal model.

We previously compared all runs with NEAMAP added and overdispersion vs catchability covariates (using maximum likelihood to calculate AIC instead of REML because we are not changing the random effects).

This work found that the models with predator attributes (mean length and number of predator species at a station) as catchability coefficients were better supported by the data than models with no catchability coefficients or with vessel effects as overdispersion parameters, similar to the findings of Ng et al. (2021), where predator mean length was selected.

Adding the predator length covariate may more directly model vessel differences in predator catch that affect stomach contents than modeling a vessel catchability covariate direclty. This was the appropach taken by Ng et al. (2021). They found that predator length covariates were strongly supported as catchability covariates (larger predators being more likely to have more prey in stomachs).

The rationale for including number of predator species is that more species “sampling” the prey field at a particular station may result in a higher encounter rate (more stations with positive bluefish prey in stomachs).

We are also interested in water temperature as a catchability covariate, because temperature affects predator feeding rate, but not all survey stations in the dataset had either surface or bottom temperature values. Rather than drop 10% of the observations (which reduced data significantly in several years), we explored different sources of surface water temperature to fill in for survey stations lacking this covariate.

Our exploration of different SST datasets found enuugh agreement with the OISST reanalysis product to use that as a first try to fill in missing data.

Here we evaluate relative fit among all covariates, as all of these models have the same input data. The SST runs were completed with this script:

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


# 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

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

MABinshore <- c(3010:3450, 3470, 3500, 3510)
GBinshore  <- c(3460, 3480, 3490, 3520:3550)

MABoffshore <- c(1010:1080, 1100:1120, 1600:1750)
GBoffshore  <- c(1090, 1130:1210, 1230, 1250)

AllEPU <- northwest_atlantic_grid %>% 
  filter(stratum_number %in% c(MAB, GB, GOM, SS)) %>% 
  select(stratum_number) %>% 
  distinct()

MABGB <- northwest_atlantic_grid %>% 
  filter(stratum_number %in% c(MAB, GB)) %>% 
  select(stratum_number) %>% 
  distinct()

MABGBinshore <- northwest_atlantic_grid %>% 
  filter(stratum_number %in% c(MABinshore, GBinshore)) %>% 
  select(stratum_number) %>% 
  distinct()

MABGBoffshore <- northwest_atlantic_grid %>% 
  filter(stratum_number %in% c(MABoffshore, GBoffshore)) %>% 
  select(stratum_number) %>% 
  distinct()

bfall <- northwest_atlantic_grid %>% 
  filter(stratum_number %in% c(bfstrata, bfoffshore)) %>% 
  select(stratum_number) %>% 
  distinct()

bfallnot <- northwest_atlantic_grid %>%
  filter(!(stratum_number %in% bfall)) %>%
  select(stratum_number) %>% 
  distinct()

bfin <- northwest_atlantic_grid %>% 
  filter(stratum_number %in% bfstrata) %>% 
  select(stratum_number) %>% 
  distinct()

bfinnot <- northwest_atlantic_grid %>%
  filter(!(stratum_number %in% bfin)) %>%
  select(stratum_number) %>% 
  distinct()

bfoff <- northwest_atlantic_grid %>% 
  filter(stratum_number %in% bfoffshore) %>% 
  select(stratum_number) %>% 
  distinct()

bfoffnot <- northwest_atlantic_grid %>%
  filter(!(stratum_number %in% bfoff)) %>%
  select(stratum_number) %>% 
  distinct()

MABGBalbinshore <- MABGBinshore %>%
  filter(!(stratum_number %in% bfstrata)) %>%
  distinct()

MABGBoffshorebigin <- MABGB %>%
  filter(!(stratum_number %in% MABGBalbinshore$stratum_number)) %>%
  distinct()

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

# Model selection options, FieldConfig default (all IID)
# Season_knots + suffix below 
# _base         No vessel overdispersion or length/number covariates  (ensure same dataset)  
# _len          Predator mean length covariate
# _no           Number of predator species covariate
# _lenno        Predator mean length and number of predator species covariates
# _eta10        Overdispersion (vessel effect) in first linear predictor (prey encounter)
# _eta11        Overdispersion (vessel effect) in both linear predictors (prey wt)

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" = AllEPU, 
                         "MABGB" = MABGB, 
                         "MABGBinshore" = MABGBinshore, 
                         "MABGBoffshore" = MABGBoffshore, 
                         "bfall" = bfall,
                         "bfallnot" = bfallnot,
                         "bfin" = bfin,
                         "bfinnot" = bfinnot,
                         "bfoff" = bfoff,
                         "bfoffnot" = bfoffnot,
                         "MABGBalbinshore" = MABGBalbinshore,
                         "MABGBoffshorebigin" = MABGBoffshorebigin))

settings = make_settings( n_x = 500, 
                          Region = "northwest_atlantic",
                          #strata.limits = list('All_areas' = 1:1e5), full area
                          strata.limits = strata.limits,
                          purpose = "index2", 
                          bias.correct = FALSE,
                          #use_anisotropy = FALSE,
                          #fine_scale = FALSE,
                          #FieldConfig = FieldConfig,
                          #RhoConfig = RhoConfig,
                          OverdispersionConfig = OverdispersionConfig
                          )


# select dataset and set directory for output

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

season <- c("fall_500_lennosst_split")

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

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


fit <- fit_model(
  settings = settings, 
  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, "/"))

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

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

season <- c("spring_500_lennosst_split")

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

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

fit = fit_model( settings = settings, 
                 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, "/"))

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

We’ll scrape the settings and parameter estimates from each folder to compare AICs:

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

# 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])
  
  
  # return model atributes 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
  )
    
    return(out)

}

# combine into one table for comparison

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

Now compare the AIC for the 500 knot models to see if including covariates or vessel effects improved the fit. This follows the model selection process outlined in Ng et al. (2021).

We have only NEAMAP data for 2020-2021 in Fall and for 2021 in Spring.

# only compare AIC for the 500 knot models
modselect.cov <- modselect %>%
  filter(n_x == 500) %>%
  #filter(str_detect(modname, "base|eta|len|_no$")) %>%
  mutate(season = case_when(str_detect(modname, "fall") ~ "Fall",
                            str_detect(modname, "spring") ~ "Spring",
                            TRUE ~ as.character(NA))) %>%
  mutate(converged2 = case_when(str_detect(converged, "no evidence") ~ "likely",
                                str_detect(converged, "is likely not") ~ "unlikely",
                                TRUE ~ as.character(NA))) %>%
  group_by(season) %>%
  mutate(deltaAIC = AIC-min(AIC)) %>%
  select(modname, season, deltaAIC, fixedcoeff,
         randomcoeff, use_anisotropy, 
         omega1, omega2, epsilon1, epsilon2, 
         beta1, beta2, AIC, converged2) %>%
  arrange(AIC)

DT::datatable(modselect.cov, rownames = FALSE, 
              options= list(pageLength = 25, scrollX = TRUE))

All models converged. Models including SST, predator mean length, and number of predator species as covariates were best supported by the data according to AIC.SST made a bigger difference in spring than in the fall model.

Below we have the results from the best fit models for each season.

Selected Model Results

Now splitting the index into many parts.

Make a lookup table:

# strata.limits <- as.list(c("AllEPU" = AllEPU, 
#                          "MABGB" = MABGB, 
#                          "MABGBinshore" = MABGBinshore, 
#                          "MABGBoffshore" = MABGBoffshore, 
#                          "bfall" = bfall,
#                          "bfallnot" = bfallnot,
#                          "bfin" = bfin,
#                          "bfinnot" = bfinnot,
#                          "bfoff" = bfoff,
#                          "bfoffnot" = bfoffnot,
#                          "MABGBalbinshore" = MABGBalbinshore,
#                          "MABGBoffshorebigin" = MABGBoffshorebigin))

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"),
                           Region  = c("AllEPU", 
                                       "MABGB", 
                                       "MABGBinshore", 
                                       "MABGBoffshore", 
                                       "bfall",
                                       "bfallnot",
                                       "bfin",
                                       "bfinnot",
                                       "bfoff",
                                       "bfoffnot",
                                       "MABGBalbinshore",
                                       "MABGBoffshorebigin"))

Fall Index

Plot individual time series

splitoutput <- read.csv("pyindex/allagg_fall_500_lennosst_split/Index.csv")

splitoutput <- splitoutput %>%
  left_join(stratlook)

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(~Region, scales = "free_y") +
  guides(colour = guide_legend(ncol=2)) +
  #theme(legend.position = c(1, 0),
   #     legend.justification = c(1, 0))
  theme(legend.position="none")

or just the indices from inshore (alb), inshore bluefish, offshore bluefish, and further out

in2off <- splitoutput %>%
  dplyr::select(Time, Region, Estimate) %>%
  tidyr::pivot_wider(names_from = Region, values_from = Estimate) %>%
  dplyr::mutate(AlbInshore = MABGBalbinshore,
                BlueInshore = bfin,
                BlueOffshore = bfoff,
                OthOffshore = MABGB - (bfoff + bfin + MABGBalbinshore),
                SumMABGB = AlbInshore + BlueInshore + BlueOffshore + OthOffshore) %>%
  dplyr::select(Time, AlbInshore, BlueInshore, BlueOffshore, OthOffshore, SumMABGB, MABGB) %>%
  tidyr::pivot_longer(!Time, names_to = "Region", values_to = "Estimate")

ggplot(in2off, aes(x=Time, y=Estimate, colour = Region)) +
  geom_point()+
  geom_line()+
  #facet_wrap(~Region) + #+ , scales = "free_y"
  #theme(legend.position = c(1, 0),
  #      legend.justification = c(1, 0))
  ggtitle("Fall Prey Index, Mid-Atlantic and Georges Bank")

or as proportions (here proportion of MABGB index).

MABGBprop <- in2off %>%
  #dplyr::filter(Region != "AllEPU") %>%
  dplyr::select(Time, Region, Estimate) %>%
  tidyr::pivot_wider(names_from = Region, values_from = Estimate) %>%
  dplyr::mutate(AlbInshoreprop = AlbInshore/MABGB,
                BlueInshoreprop = BlueInshore/MABGB,
                BlueOffshoreprop = BlueOffshore/MABGB,
                OthOffshoreprop = OthOffshore/MABGB) %>%
  tidyr::pivot_longer(!Time, names_to = "Region", values_to = "Estimate") %>%
  dplyr::filter(Region %in% c("AlbInshoreprop", "BlueInshoreprop", "BlueOffshoreprop",
                              "OthOffshoreprop"))
  

ggplot(MABGBprop, aes(x=Time, y=Estimate, colour = Region)) +
  geom_point()+
  geom_line()+
  #facet_wrap(~Region) + #+ , scales = "free_y"
  #theme(legend.position = c(1, 0),
  #      legend.justification = c(1, 0))
  ggtitle("Fall Prey Index as proportion of Mid-Atlantic and Georges Bank")

Fall predicted ln-density

Fall density maps with covariates

Spring Index

Plot individual time series

splitoutput <- read.csv("pyindex/allagg_spring_500_lennosst_split/Index.csv")

splitoutput <- splitoutput %>%
  left_join(stratlook)

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(~Region, scales = "free_y") +
  guides(colour = guide_legend(ncol=2)) +
  #theme(legend.position = c(1, 0),
   #     legend.justification = c(1, 0))
  theme(legend.position="none")

or just the indices from inshore (alb), inshore bluefish, offshore bluefish, and further out

in2off <- splitoutput %>%
  dplyr::select(Time, Region, Estimate) %>%
  tidyr::pivot_wider(names_from = Region, values_from = Estimate) %>%
  dplyr::mutate(AlbInshore = MABGBalbinshore,
                BlueInshore = bfin,
                BlueOffshore = bfoff,
                OthOffshore = MABGB - (bfoff + bfin + MABGBalbinshore),
                SumMABGB = AlbInshore + BlueInshore + BlueOffshore + OthOffshore) %>%
  dplyr::select(Time, AlbInshore, BlueInshore, BlueOffshore, OthOffshore, SumMABGB, MABGB) %>%
  tidyr::pivot_longer(!Time, names_to = "Region", values_to = "Estimate")

ggplot(in2off, aes(x=Time, y=Estimate, colour = Region)) +
  geom_point()+
  geom_line()+
  #facet_wrap(~Region) + #+ , scales = "free_y"
  #theme(legend.position = c(1, 0),
  #      legend.justification = c(1, 0))
  ggtitle("Spring Prey Index, Mid-Atlantic and Georges Bank")

or as proportions (here proportion of MABGB index).

MABGBprop <- in2off %>%
  #dplyr::filter(Region != "AllEPU") %>%
  dplyr::select(Time, Region, Estimate) %>%
  tidyr::pivot_wider(names_from = Region, values_from = Estimate) %>%
  dplyr::mutate(AlbInshoreprop = AlbInshore/MABGB,
                BlueInshoreprop = BlueInshore/MABGB,
                BlueOffshoreprop = BlueOffshore/MABGB,
                OthOffshoreprop = OthOffshore/MABGB) %>%
  tidyr::pivot_longer(!Time, names_to = "Region", values_to = "Estimate") %>%
  dplyr::filter(Region %in% c("AlbInshoreprop", "BlueInshoreprop", "BlueOffshoreprop",
                              "OthOffshoreprop"))
  

ggplot(MABGBprop, aes(x=Time, y=Estimate, colour = Region)) +
  geom_point()+
  geom_line()+
  #facet_wrap(~Region) + #+ , scales = "free_y"
  #theme(legend.position = c(1, 0),
  #      legend.justification = c(1, 0))
  ggtitle("Spring Prey Index as proportion of Mid-Atlantic and Georges Bank")

Spring predicted ln-density

Spring density maps with covariates

All results in the respective pyindex/allagg_fall_500_lennosst_split and allagg_spring_500_lennsst_split folders.

The full results are on google drive rather than github to save space.

Still to do:

  • index within 3 miles of shore and outside that (highest priority)

  • investigate other SST filling sources if time

    • CTD casts from survey not meeting criteria for proximity to station
    • underway measurements on Bigelow and possibly Albatross
    • AVHRR satellite data

References

Ng, E. L., Deroba, J. J., Essington, T. E., Grüss, A., Smith, B. E., and Thorson, J. T. 2021. Predator stomach contents can provide accurate indices of prey biomass. ICES Journal of Marine Science, 78: 1146–1159. https://doi.org/10.1093/icesjms/fsab026 (Accessed 1 September 2021).