Model comparisons led us to the best model fit using mean predator length, number of predator species, and SST at a survey station as catchability covariates.

The last step is to define spatial partitions that reflect several definitions of “inshore” and “offshore” for the stock assessment inputs.

First we define a partition that is the MAB and GB areas only as the GOM is not relevant to the bluefish assessment (yet). This is called MABGB. Within this partition,

  1. Survey inshore vs offshore to evaluate availability to the survey index. Strata partitions include:
    • Albatross inshore stations
    • Bigelow inshore bluefish index stations
    • offshore bluefish index stations (added this year)
    • offshore non-bluefish stations
  2. Recreational fishery inshore vs offshore to evaluate availability to the MRIP CPUE index. Strata partitions include
    • shoreline to 3 miles out
    • offshore of 3 miles

Survey strata definitions are built into VAST already.

The area within and outside 3 miles of shore within MABGB needs to be defined. At the moment, 3 nm is approximated as 5.556 km

## Northeast:

# set bounding boxes
neus.xmin=-77
neus.xmax=-65
neus.ymin=35
neus.ymax=45

# high resolution coastline
usamap <- rnaturalearth::ne_countries(scale = "large", country = "united states of america", returnclass = "sf")[1] %>% 
  sf::st_cast("MULTILINESTRING") # get basic map of the country 

neus.bbox1 <- sf::st_set_crs(sf::st_as_sf(as(raster::extent(neus.xmin, neus.xmax, neus.ymin, neus.ymax), "SpatialPolygons")), sf::st_crs(usamap))
neus.bbox2 <- sf::st_set_crs(sf::st_as_sf(as(raster::extent(-78, -74, 42, 45), "SpatialPolygons")), sf::st_crs(usamap)) # smaller bounding box to get rid of extra lines on the map 

# just the NEUS coastline high res

neuscoast <- usamap %>% 
  sf::st_intersection(neus.bbox1) %>%  
  sf::st_difference(neus.bbox2) # gets rid of extra non coastal line 

#plot(neuscoast)

# add a 5.556 km (3 nautical mi) buffer around coastline

neuscoast_buff_3nm  <-  sf::st_buffer(neuscoast, dist = 5556)

#plot(neuscoast_buff_3nm)

# intersect buffer with the current FishStatsUtils::northwest_atlantic_grid
# make northwest atlantic grid into sf object
nwagrid_sf  <-  sf::st_as_sf(FishStatsUtils::northwest_atlantic_grid, coords = c("Lon","Lat")) %>%
  sf::st_set_crs(sf::st_crs(neuscoast))

# intersect, rearrange in same format as nwatl grid, and save
coast3nmbuff <- sf::st_intersection(nwagrid_sf,neuscoast_buff_3nm) %>%
  dplyr::mutate(Lon = as.numeric(sf::st_coordinates(.)[,1]),
                Lat = as.numeric(sf::st_coordinates(.)[,2])) %>%
  sf::st_set_geometry(NULL) %>%
  dplyr::select(-featurecla) %>%
  dplyr::select(stratum_number, Lon, Lat, everything())

write_rds(coast3nmbuff, here("spatialdat","neus_coast3nmbuff.rds"))

fedwaters <- setdiff(FishStatsUtils::northwest_atlantic_grid, coast3nmbuff)

Plot the 3 mile grid locations

ggplot(data = ecodata::coast) +
  geom_sf() + 
  geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat), size=0.05, alpha=0.1) +
  geom_point(data = coast3nmbuff, aes(x = Lon, y = Lat), size=0.05, colour = "green",  alpha=0.1) +
  coord_sf(xlim = c(-78, -65.5), ylim = c(35, 45))+
  ggtitle("State waters (3 nm)")

Plot the fed waters outside 3 miles

ggplot(data = ecodata::coast) +
  geom_sf() + 
  geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat), size=0.05, alpha=0.1) +
  geom_point(data = fedwaters, aes(x = Lon, y = Lat), size=0.05, colour = "green",  alpha=0.1) +
  coord_sf(xlim = c(-78, -65.5), ylim = c(35, 45))+
  ggtitle("Federal waters")

Zoom in

ggplot(data = ecodata::coast) +
  geom_sf() + 
  geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat), size=0.05, alpha=0.1) +
  geom_point(data = coast3nmbuff, aes(x = Lon, y = Lat), size=0.05, colour = "green",  alpha=0.1) +
  #coord_sf(xlim = c(-78, -65.5), ylim = c(35, 45))
  coord_sf(xlim = c(-74.5, -69), ylim = c(38.5, 42.5)) +
  ggtitle("Zoomed State waters (3 nm)")

ggplot(data = ecodata::coast) +
  geom_sf() + 
  geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat), size=0.05, alpha=0.1) +
  geom_point(data = fedwaters, aes(x = Lon, y = Lat), size=0.05, colour = "green",  alpha=0.1) +
  #coord_sf(xlim = c(-78, -65.5), ylim = c(35, 45))
  coord_sf(xlim = c(-74.5, -69), ylim = c(38.5, 42.5))+
  ggtitle("Zoomed federal waters")

Definitions for strata:

#current bluefish assessment strata are all Bigelow inshore strata MAB-GB
bfinshore <- c(3020, 3050, 3080, 3110, 3140, 3170, 3200, 3230, 
              3260, 3290, 3320, 3350, 3380, 3410, 3440, 3450, 3460)

bfinshoregrid <-  FishStatsUtils::northwest_atlantic_grid %>%
  filter(stratum_number %in% bfinshore)
  
  
# from Tony's 8 March presentation, minus the inshore in CCBay
bfoffshore <- c(1010, 1730, 1690, 1650, 1050, 1060, 1090, 1100, 1250, 1200, 1190, 1610)

bfoffshoregrid <-  FishStatsUtils::northwest_atlantic_grid %>%
  filter(stratum_number %in% bfoffshore)

#from mskeyrun vignette, EPU based on survey strata, replace built in VAST EPU
#https://noaa-edab.github.io/ms-keyrun/articles/GBSurveySet.html

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)

MABGBgrid <-  FishStatsUtils::northwest_atlantic_grid %>%
  filter(stratum_number %in% c(MAB, GB))

albinshoregrid <- MABGBgrid %>%
  filter(stratum_number>2999 & stratum_number<3999) %>% #inshore
  anti_join(bfinshoregrid)

othoffshoregrid <- MABGBgrid %>%
  anti_join(bind_rows(albinshoregrid, bfinshoregrid, bfoffshoregrid))

statewatersgrid <- coast3nmbuff %>%
  inner_join(MABGBgrid)

fedwatersgrid <- fedwaters %>%
  inner_join(MABGBgrid)

Plot areas

ggplot(data = ecodata::coast) +
  geom_sf() + 
  geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat), size=0.05, alpha=0.1) +
  geom_point(data = albinshoregrid, aes(x = Lon, y = Lat), size=0.05, colour = "green",  alpha=0.1) +
  coord_sf(xlim = c(-78, -65.5), ylim = c(35, 45)) +
  ggtitle("Albatross inshore strata")

ggplot(data = ecodata::coast) +
  geom_sf() + 
  geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat), size=0.05, alpha=0.1) +
  geom_point(data = bfinshoregrid, aes(x = Lon, y = Lat), size=0.05, colour = "green",  alpha=0.1) +
  coord_sf(xlim = c(-78, -65.5), ylim = c(35, 45)) +
  ggtitle("Bluefish inshore survey strata")

ggplot(data = ecodata::coast) +
  geom_sf() + 
  geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat), size=0.05, alpha=0.1) +
  geom_point(data = bfoffshoregrid, aes(x = Lon, y = Lat), size=0.05, colour = "green",  alpha=0.1) +
  coord_sf(xlim = c(-78, -65.5), ylim = c(35, 45)) +
  ggtitle("Bluefish offshore survey strata")

ggplot(data = ecodata::coast) +
  geom_sf() + 
  geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat), size=0.05, alpha=0.1) +
  geom_point(data = othoffshoregrid, aes(x = Lon, y = Lat), size=0.05, colour = "green",  alpha=0.1) +
  coord_sf(xlim = c(-78, -65.5), ylim = c(35, 45)) +
  ggtitle("Other offshore strata")

ggplot(data = ecodata::coast) +
  geom_sf() + 
  geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat), size=0.05, alpha=0.1) +
  geom_point(data =statewatersgrid, aes(x = Lon, y = Lat), size=0.05, colour = "green",  alpha=0.1) +
  coord_sf(xlim = c(-78, -65.5), ylim = c(35, 45)) +
  ggtitle("State waters")

ggplot(data = ecodata::coast) +
  geom_sf() + 
  geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat), size=0.05, alpha=0.1) +
  geom_point(data = fedwatersgrid, aes(x = Lon, y = Lat), size=0.05, colour = "green",  alpha=0.1) +
  coord_sf(xlim = c(-78, -65.5), ylim = c(35, 45)) +
  ggtitle("Federal waters")

Run 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

#current bluefish assessment strata are all Bigelow inshore strata MAB-GB
bfinshore <- c(3020, 3050, 3080, 3110, 3140, 3170, 3200, 3230, 
               3260, 3290, 3320, 3350, 3380, 3410, 3440, 3450, 3460)

bfinshoregrid <-  FishStatsUtils::northwest_atlantic_grid %>%
  filter(stratum_number %in% bfinshore)

# from Tony's 8 March presentation, minus the inshore in CCBay
bfoffshore <- c(1010, 1730, 1690, 1650, 1050, 1060, 1090, 1100, 1250, 1200, 1190, 1610)

bfoffshoregrid <-  FishStatsUtils::northwest_atlantic_grid %>%
  filter(stratum_number %in% bfoffshore)

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


#from mskeyrun vignette, EPU based on survey strata, replace built in VAST EPU
#https://noaa-edab.github.io/ms-keyrun/articles/GBSurveySet.html

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)

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

MABGBgrid <-  FishStatsUtils::northwest_atlantic_grid %>%
  filter(stratum_number %in% c(MAB, GB))

MABGB <- MABGBgrid %>% 
  select(stratum_number) %>% 
  distinct()

albinshoregrid <- MABGBgrid %>%
  filter(stratum_number>2999 & stratum_number<3999) %>% #inshore
  anti_join(bfinshoregrid)

albinshore <- albinshoregrid %>%
  select(stratum_number) %>% 
  distinct() 

othoffshoregrid <- MABGBgrid %>%
  anti_join(bind_rows(albinshoregrid, bfinshoregrid, bfoffshoregrid))

othoffshore <- othoffshoregrid %>%
  select(stratum_number) %>% 
  distinct() 

coast3nmbuff <- readRDS(here("spatialdat/neus_coast3nmbuff.rds"))

fedwaters <- setdiff(FishStatsUtils::northwest_atlantic_grid, coast3nmbuff)

statewatersgrid <- coast3nmbuff %>%
  inner_join(MABGBgrid)

fedwatersgrid <- fedwaters %>%
  inner_join(MABGBgrid)


# 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, 
                         "albinshore" = albinshore,
                         "bfinshore" = bfinshore,
                         "bfoffshore" = bfoffshore,                         
                         "bfall" = bfall,
                         "othoffshore" = othoffshore
                         ))

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_split1")

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_split1")

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, "/"))

To actually use the state vs federal waters split, I need to make a user defined extrapolation grid that has these divisions as strata. So the model will be run twice, once with the standard grid and strata-based divisions, and once with the “new” grid that defines state vs federal waters as strata.

Make the new grid:

coast3nmbuff <- readRDS(here("spatialdat/neus_coast3nmbuff.rds"))

coast3nmbuffst <- coast3nmbuff %>%
  dplyr::mutate(strat2 = 1) %>% #state waters = 1
  dplyr::right_join(FishStatsUtils::northwest_atlantic_grid) %>%
  dplyr::mutate(strat2 = replace_na(strat2, 2)) %>% #replace NA with 2 for fed waters
  dplyr::mutate(strat2 = case_when(!stratum_number %in% c(MAB, GB) ~ 0, #ignore outside MABGB
                                   TRUE ~ strat2)) %>%
  dplyr::mutate(stratum_number = strat2) %>%
  dplyr::select(-strat2)
  
saveRDS(coast3nmbuffst, file = here("spatialdat/user_region.rds"))

This is the script for the second run with the new grid.

# 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

user_region <- readRDS(here::here("spatialdat/user_region.rds"))
names(user_region)[4] <- "Area_km2"

MABGBstate <- 1
MABGBfed <- 2
MABGB <- c(1,2)

# 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 <- data.frame('STRATA' = c("MABGBstate", "MABGBfed"),
                  "stratum_number" = c(1, 2))
                            

settings = make_settings( n_x = 500, 
                          Region = "User",
                          #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_3mi")

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,
  input_grid=user_region,
  working_dir = paste0(working_dir, "/"))

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

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

season <- c("spring_500_lennosst_3mi")

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,
                 input_grid=user_region,
                 working_dir = paste0(working_dir, "/"))

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

Selected Model Results: 3 mi split

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"),
                           Region  = c("state", 
                                       "federal"))

Fall Index

Plot individual time series

splitoutput <- read.csv("pyindex/allagg_fall_500_lennosst_3mi/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")

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_3mi/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:

  • fix 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