VAST aggregate prey model testing different strata.limits
NEFSC and NEAMAP surveys, new predator list, all covariates
Sarah Gaichas
15 July, 2022
Investigating the stability of index estimates using different strata limits definitions. Ideally, the index for our areas of interest will be robust to the selection of overall model strata.limits (within reason).
The largest possible spatial domain is the full Northwest Atlantic grid in VAST, which is the full set of NEFSC survey strata including those not routinely sampled anymore south of Cape Hatteras. The smallest possible spatial domain for the bluefish assessment is the MABGB survey strata. I would rather not model at this level if we want to use the same dataset and VAST model to derive an index for the Gulf of Maine in addition to Mid Atlantic and Georges Bank for State of the Ecosystem reporting. But we’ll look at both extremes as the overall strata.limits here and compare index results.
Why am I opening this can of worms: the fall model converges using the full set of “EPU” strata in the VAST northwest_atlantic_grid, and when using all survey strata including those outside the defined EPUs, but does not converge when limited to the survey strata defined AllEPU strata limits, which differ only a little spatially. Spring converges in both cases. Any change in the model spatial domain overall will result in different parameter estimates. But it is the index we are interested in.
Model comparisons evaluating fit with different catchability coefficents, etc, for the full model converge in fall because we are using the entire survey grid that extends beyond either EPU definition (South of Cape Hatteras and maybe also a bit north), and then splitting to our areas of interest. As noted during analysis of the bias corrected run,
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. So I was incorrect in thinking I was using the same strata.limits, and was correct in my understanding of how the model works. I just failed to notice I was actually comparing the full range of survey strata with the truncated range in AllEPU as the model outer limits.
So the question is whether we get a similar index using the entire survey grid and partitioning results to AllEPU or some other strata, in which case we go ahead and run on the full survey grid and don’t worry about it.
Here we compare indices estimated for each level of
strata.limits defined in the VAST run.
Definitions for strata used in the script
bluefishdiet/VASTunivariate_bfp_teststratalimits.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()
MABGBGOM <- northwest_atlantic_grid %>%
filter(stratum_number %in% c(MAB, GB, GOM)) %>%
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()Visuals of what the full model domain looks like with strata.limits
defined using all of these groupings that extend to all survey strata
(model name ends in _split1):
split1 spatial domain
And this is the model domain when strata.limits are only survey
strata-defined EPUs (model name ends in _AllEPU):
AllEPU spatial domain
And finally just for testing, the model domain when strata.limits
include only the MAB and GB survey strata-defined EPUs (model name ends
in _MABGB):
MABGB spatial domain
We are now back to bias correction turned off, and instead running this script for a single set of strata.limits to compare indices generated:
# 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()
MABGBGOM <- northwest_atlantic_grid %>%
filter(stratum_number %in% c(MAB, GB, GOM)) %>%
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(MABGBGOM
#"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_MABGBGOM")
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_MABGBGOM")
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, comparing strata limit definitions
Compare full model (all survey strata including below Hatteras) partitioned results to those estimated only from AllEPU, MABGB:
moddir <- c("pyindex/allagg_fall_500_lennosst_split1",
"pyindex/allagg_fall_500_lennosst_AllEPU",
"pyindex/allagg_fall_500_lennosst_MABGB",
"pyindex/allagg_spring_500_lennosst_split1",
"pyindex/allagg_spring_500_lennosst_AllEPU",
"pyindex/allagg_spring_500_lennosst_MABGB")
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"))
# function to apply extracting info
getmodinfo.index <- function(d.name){
# read settings
modpath <- stringr::str_split(d.name, "/", simplify = TRUE)
modname <- modpath[length(modpath)]
settings <- read.table(file.path(d.name, "settings.txt"), comment.char = "",
fill = TRUE, header = FALSE)
n_x <- as.numeric(as.character(settings[(which(settings[,1]=="$n_x")+1),2]))
grid_size_km <- as.numeric(as.character(settings[(which(settings[,1]=="$grid_size_km")+1),2]))
max_cells <- as.numeric(as.character(settings[(which(settings[,1]=="$max_cells")+1),2]))
use_anisotropy <- as.character(settings[(which(settings[,1]=="$use_anisotropy")+1),2])
fine_scale <- as.character(settings[(which(settings[,1]=="$fine_scale")+1),2])
bias.correct <- as.character(settings[(which(settings[,1]=="$bias.correct")+1),2])
#FieldConfig
if(settings[(which(settings[,1]=="$FieldConfig")+1),1]=="Component_1"){
omega1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+2),2])
omega2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+3),1])
epsilon1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+4),2])
epsilon2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+5),1])
beta1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+6),2])
beta2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+7),1])
}
if(settings[(which(settings[,1]=="$FieldConfig")+1),1]=="Omega1"){
omega1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+3),1])
omega2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+4),1])
epsilon1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+3),2])
epsilon2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+4),2])
beta1 <- "IID"
beta2 <- "IID"
}
#RhoConfig
rho_beta1 <- as.numeric(as.character(settings[(which(settings[,1]=="$RhoConfig")+3),1]))
rho_beta2 <- as.numeric(as.character(settings[(which(settings[,1]=="$RhoConfig")+3),2]))
rho_epsilon1 <- as.numeric(as.character(settings[(which(settings[,1]=="$RhoConfig")+4),1]))
rho_epsilon2 <- as.numeric(as.character(settings[(which(settings[,1]=="$RhoConfig")+4),2]))
# read parameter estimates, object is called parameter_Estimates
load(file.path(d.name, "parameter_estimates.RData"))
AIC <- parameter_estimates$AIC[1]
converged <- parameter_estimates$Convergence_check[1]
fixedcoeff <- unname(parameter_estimates$number_of_coefficients[2])
randomcoeff <- unname(parameter_estimates$number_of_coefficients[3])
index <- read.csv(file.path(d.name, "Index.csv"))
# return model attributes as a dataframe
out <- data.frame(modname = modname,
n_x = n_x,
grid_size_km = grid_size_km,
max_cells = max_cells,
use_anisotropy = use_anisotropy,
fine_scale = fine_scale,
bias.correct = bias.correct,
omega1 = omega1,
omega2 = omega2,
epsilon1 = epsilon1,
epsilon2 = epsilon2,
beta1 = beta1,
beta2 = beta2,
rho_epsilon1 = rho_epsilon1,
rho_epsilon2 = rho_epsilon2,
rho_beta1 = rho_beta1,
rho_beta2 = rho_beta2,
AIC = AIC,
converged = converged,
fixedcoeff = fixedcoeff,
randomcoeff = randomcoeff,
index = index
)
return(out)
}
modcompare <- purrr::map_dfr(moddir, getmodinfo.index)Fall results
Which models converged?
fallconverged <- modcompare %>%
dplyr::filter(grepl("fall", modname)) %>%
dplyr::select(modname, converged) %>%
dplyr::distinct()
knitr::kable(fallconverged)| modname | converged |
|---|---|
| allagg_fall_500_lennosst_split1 | There is no evidence that the model is not converged |
| allagg_fall_500_lennosst_AllEPU | The model is likely not converged |
| allagg_fall_500_lennosst_MABGB | There is no evidence that the model is not converged |
Fall Index
Plot index time series with standard errors by model:
falloutput <- modcompare %>%
dplyr::filter(grepl("fall", modname)) %>%
dplyr::left_join(stratlook, by = c("index.Stratum" = "Stratum")) %>%
dplyr::mutate(Region = case_when(grepl("MABGB", modname) ~ "MABGB",
TRUE ~ Region)) %>%
dplyr::filter(Region %in% c("AllEPU", "MABGB"))
ggplot(falloutput, aes(x=index.Time, y=index.Estimate, colour=modname)) +
geom_errorbar(aes(ymin=index.Estimate+index.Std..Error.for.Estimate, ymax=index.Estimate-index.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="bottom")
Spring results
springconverged <- modcompare %>%
dplyr::filter(grepl("spring", modname)) %>%
dplyr::select(modname, converged) %>%
dplyr::distinct()
knitr::kable(springconverged)| modname | converged |
|---|---|
| allagg_spring_500_lennosst_split1 | There is no evidence that the model is not converged |
| allagg_spring_500_lennosst_AllEPU | There is no evidence that the model is not converged |
| allagg_spring_500_lennosst_MABGB | There is no evidence that the model is not converged |
Spring Index
Plot index time series with standard errors by model:
springoutput <- modcompare %>%
dplyr::filter(grepl("spring", modname)) %>%
dplyr::left_join(stratlook, by = c("index.Stratum" = "Stratum")) %>%
dplyr::mutate(Region = case_when(grepl("MABGB", modname) ~ "MABGB",
TRUE ~ Region)) %>%
dplyr::filter(Region %in% c("AllEPU", "MABGB"))
ggplot(springoutput, aes(x=index.Time, y=index.Estimate, colour=modname)) +
geom_errorbar(aes(ymin=index.Estimate+index.Std..Error.for.Estimate, ymax=index.Estimate-index.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="bottom")
OK I feel better about this. All indices look very close, and well
within standard error ranges. Perhaps best to continue with the full set
of survey strata (here ..._split1 models), since those
models converge consistently. It seems reasonable that strata limits at
the MABGB level would depart a bit more from the full model narrowed
down to MABGB as they are excluding data from the full spatial
domain.
Last thing, what difference did bias correction make on the index? The above are all not bias corrected because that takes four times as long to run. However we can compare the split1 models to the bias corrected models because they were run on exactly the same footprint. Here is the difference:
Bias corrected vs not
moddir2 <- c("pyindex/allagg_fall_500_lennosst_split1",
"pyindex/allagg_fall_500_lennosst_split_biascorrect",
"pyindex/allagg_spring_500_lennosst_split1",
"pyindex/allagg_spring_500_lennosst_split_biascorrect")
modcompare2 <- purrr::map_dfr(moddir2, getmodinfo.index)Fall results
Which models converged?
fallconverged2 <- modcompare2 %>%
dplyr::filter(grepl("fall", modname)) %>%
dplyr::select(modname, converged) %>%
dplyr::distinct()
knitr::kable(fallconverged2)| modname | converged |
|---|---|
| allagg_fall_500_lennosst_split1 | There is no evidence that the model is not converged |
| allagg_fall_500_lennosst_split_biascorrect | There is no evidence that the model is not converged |
Fall Index
Plot index time series with standard errors by model:
falloutput2 <- modcompare2 %>%
dplyr::filter(grepl("fall", modname)) %>%
dplyr::left_join(stratlook, by = c("index.Stratum" = "Stratum")) %>%
#dplyr::mutate(Region = case_when(grepl("MABGB", modname) ~ "MABGB",
# TRUE ~ Region)) %>%
dplyr::filter(Region %in% c("AllEPU", "MABGB",
"bfin", "bfoff",
"MABGBalbinshore"))
ggplot(falloutput2, aes(x=index.Time, y=index.Estimate, colour=modname)) +
geom_errorbar(aes(ymin=index.Estimate+index.Std..Error.for.Estimate, ymax=index.Estimate-index.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="bottom")
Spring results
Which models converged?
springconverged2 <- modcompare2 %>%
dplyr::filter(grepl("spring", modname)) %>%
dplyr::select(modname, converged) %>%
dplyr::distinct()
knitr::kable(springconverged2)| modname | converged |
|---|---|
| allagg_spring_500_lennosst_split1 | There is no evidence that the model is not converged |
| allagg_spring_500_lennosst_split_biascorrect | There is no evidence that the model is not converged |
Spring index
Plot index time series with standard errors by model:
springoutput2 <- modcompare2 %>%
dplyr::filter(grepl("spring", modname)) %>%
dplyr::left_join(stratlook, by = c("index.Stratum" = "Stratum")) %>%
#dplyr::mutate(Region = case_when(grepl("MABGB", modname) ~ "MABGB",
# TRUE ~ Region)) %>%
dplyr::filter(Region %in% c("AllEPU", "MABGB",
"bfin", "bfoff",
"MABGBalbinshore"))
ggplot(springoutput2, aes(x=index.Time, y=index.Estimate, colour=modname)) +
geom_errorbar(aes(ymin=index.Estimate+index.Std..Error.for.Estimate, ymax=index.Estimate-index.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="bottom")
Bias correction tends to inflate the index a bit.
This probably doesn’t matter if it is used as a covariate in an assessment model.