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 model has been partitioned into 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, while the full model is AllEPU. 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 results presented here pertain only to survey strata; I’ll consider the 3 mile split separately in a different model run. I can’t find a way to have VAST do all the partitions at once, because the 3 mile split cuts across survey stratum boundaries.

Definitions for strata used in the script bluefishdiet/VASTunivariate_bfp_allsurvs_lencovSST_inoffsplit.R:

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

Attempts to rewrite these strata to a subset of only the needed ones failed! The model did not converge! I need to better understand how stratum definitions affect model estimation. I thought–apparently incorrectly–that they didn’t. This configuration, where each strata set also has a complement included, results in model convergence.

Bias correction is based on (Thorson and Kristensen, 2016).

Run this script, the same used to test catchability covariates, but now has all three covariates (predator mean length, number of predator species, and SSTfill) included and bias.correct=TRUE:

# 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 = TRUE,
                          #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_biascorrect")

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

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

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

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

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

Model Results, Bias-corrected survey strata

Fall results

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 (bias corrected) with standard errors:

splitoutput <- read.csv("pyindex/allagg_fall_500_lennosst_split_biascorrect/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, Std..Error.for.Estimate) %>%
  tidyr::pivot_longer(c(Estimate, Std..Error.for.Estimate), names_to = "Var") %>%
  dplyr::group_by(Var) %>%
  tidyr::pivot_wider(names_from = Region, values_from = value) %>%
  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(!c(Time, Var), names_to = "Region", values_to = "value") %>%
  tidyr::pivot_wider(names_from = "Var", values_from = "value")

ggplot(in2off, 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"
  #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

Fall Diagnostics

diagplots <- c("Data_and_knots",
               "Data_by_year",
               "quantile_residuals",
               "quantile_residuals_on_map",
               "Aniso",
               "center_of_gravity")

for(p in diagplots){
  
    cat("  \n####",  as.character(p),"  \n")
    cat(paste0("![](pyindex/allagg_fall_500_lennosst_split_biascorrect/",
                      p,
                      ".png)")) 
    cat("  \n")   
    
  }

Data_and_knots

Data_by_year

quantile_residuals

quantile_residuals_on_map

Aniso

center_of_gravity

Spring results

Spring Index

Plot individual time series

splitoutput <- read.csv("pyindex/allagg_spring_500_lennosst_split_biascorrect/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, Std..Error.for.Estimate) %>%
  tidyr::pivot_longer(c(Estimate, Std..Error.for.Estimate), names_to = "Var") %>%
  dplyr::group_by(Var) %>%
  tidyr::pivot_wider(names_from = Region, values_from = value) %>%
  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(!c(Time, Var), names_to = "Region", values_to = "value") %>%
  tidyr::pivot_wider(names_from = "Var", values_from = "value")

ggplot(in2off, 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"
  #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

Spring Diagnostics

diagplots <- c("Data_and_knots",
               "Data_by_year",
               "quantile_residuals",
               "quantile_residuals_on_map",
               "Aniso",
               "center_of_gravity")

for(p in diagplots){
  
    cat("  \n####",  as.character(p),"  \n")
    cat(paste0("![](pyindex/allagg_spring_500_lennosst_split_biascorrect/",
                      p,
                      ".png)")) 
    cat("  \n")   
    
  }

Data_and_knots

Data_by_year

quantile_residuals

quantile_residuals_on_map

Aniso

center_of_gravity

All results in the respective pyindex/allagg_fall_500_lennosst_split_biascorrect and allagg_spring_500_lennsst_splitbiascorrect 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)

References

Thorson, J.T., and Kristensen, K.
2016. Implementing a generic method for bias correction in statistical models using random effects, with spatial and population dynamics examples. Fisheries Research 175:66–74. https://doi.org/10.1016/j.fishres.2015.11.016.