ZoopCOG
Sarah Gaichas
2024-10-30
Pull COG from “best” zooplankton models
See Copepod Model Results for the model stats, time series, and density outputs
These models are for all manner of copepods, plus zooplankton volume and Euphausiids
Functions to get the COG out and plot it with geom_gls trend assessment
# put this in a package for petes sake
getcogVAST <- function(d.name){
fit <- VAST::reload_model(readRDS(paste0(d.name,"/fit.rds")))
dir.create(paste0(d.name,"/test"))
cogout <- FishStatsUtils::plot_range_index(Sdreport = fit$parameter_estimates$SD,
Report = fit$Report,
TmbData = fit$data_list,
year_labels = as.numeric(fit$year_labels),
years_to_plot = fit$years_to_plot,
Znames = colnames(fit$data_list$Z_xm),
PlotDir = paste0(d.name,"/test")) #already have plots, will delete this directory
saveRDS(cogout, paste0(d.name,"/cogout.rds"))
unlink(paste0(d.name,"/test"), recursive=TRUE) #removes directory with unneeded plots
}
plotcog <- function(d.name){
cogout <- readRDS(paste0(d.name,"/cogout.rds"))
modpath <- unlist(str_split(d.name, pattern = "/"))
modname <- modpath[length(modpath)]
cogdat <- as.data.frame(cogout$COG_Table) |>
dplyr::mutate(direction = ifelse(m==1, "Eastward", "Northward"))
ggplot2::ggplot(cogdat, ggplot2::aes(x = Year, y = COG_hat)) +
ggplot2::geom_point() +
ecodata::geom_gls() +
ggplot2::labs(y = "Center of gravity, km") +
ggplot2::facet_wrap(~direction, scales = "free_y") +
ecodata::theme_facet()+
ggplot2::ggtitle(modname)
}Functions to get the “best” model directories and apply the function to them
modtable <- function(moddirs){
# apply getmodinfo function to inout directories
modcompare <- purrr::map_dfr(moddirs, getmodinfo)
modselect <- modcompare |>
dplyr::mutate(season = dplyr::case_when(stringr::str_detect(modname, "_fall_") ~ "Fall",
stringr::str_detect(modname, "spring") ~ "Spring",
stringr::str_detect(modname, "_all_") ~ "Annual",
TRUE ~ as.character(NA))) |>
dplyr::mutate(converged2 = dplyr::case_when(stringr::str_detect(converged, "no evidence") ~ "likely",
stringr::str_detect(converged, "is likely not") ~ "unlikely",
TRUE ~ as.character(NA))) |>
dplyr::mutate(copegroup = stringr::str_extract(modname, "[^_]+")) |>
#dplyr::mutate(modname = str_extract(modname, '(?<=allagg_).*')) |>
dplyr::group_by(copegroup, season) |>
dplyr::mutate(deltaAIC = AIC-min(AIC)) |>
dplyr::select(copegroup, modname, season, deltaAIC, fixedcoeff,
randomcoeff, use_anisotropy,
omega1, omega2, epsilon1, epsilon2,
beta1, beta2, AIC, converged2) |>
dplyr::arrange(copegroup, season, AIC)
return(modselect)
}
# 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
if(file.exists(file.path(d.name, "parameter_estimates.RData"))) {
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])
}else if(file.exists(file.path(d.name, "parameter_estimates.txt"))){
reptext <- readLines(file.path(d.name, "parameter_estimates.txt"))
AIC <- as.double(reptext[grep(reptext, pattern = "AIC")+1])
converged <- reptext[grep(reptext, pattern = "Convergence_check")+1]
fixedcoeff <- as.integer(stringr::str_split(reptext[grep(reptext, pattern = "number_of_coefficients")+2],
boundary("word"))[[1]][2])
randomcoeff <- as.integer(stringr::str_split(reptext[grep(reptext, pattern = "number_of_coefficients")+2],
boundary("word"))[[1]][3])
}else{
AIC <- NA_real_
converged <- NA_character_
fixedcoeff <- NA_integer_
randomcoeff <- NA_integer_
}
#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)
}
# function to apply extracting info
getmodindex <- function(d.name){
# read settings
modpath <- stringr::str_split(d.name, "/", simplify = TRUE)
modname <- modpath[length(modpath)]
if(file.exists(file.path(d.name,"Index.csv"))){
index <- read.csv(file.path(d.name, "Index.csv"))
}else{
stopifnot()
}
# return model indices as a dataframe
out <- data.frame(modname = modname,
index
)
return(out)
}Apply all that
# from each output folder in pyindex,
outdir <- here::here("pyindex")
moddirs <- list.dirs(outdir)
moddirs <- moddirs[-1]
modselect <- modtable(moddirs)
# lets only look at indices for converged models
modcompare_conv <- modselect |>
dplyr::ungroup() |>
dplyr::filter(converged2 == "likely") |>
dplyr::select(modname) |>
as.vector() |>
unname() |>
unlist()
# or only converged and best AIC
modcompare_best <- modselect |>
dplyr::ungroup() |>
dplyr::filter(converged2 == "likely") |>
dplyr::distinct(copegroup, season, .keep_all = T) |> # assumes orderd by dAIC!
dplyr::select(modname) |>
as.vector() |>
unname() |>
unlist()
moddirs_conv <- moddirs[grepl(sprintf("\\.*(%s)$", paste(modcompare_conv, collapse = '|')), moddirs)]
moddirs_best <- moddirs[grepl(sprintf("\\.*(%s)$", paste(modcompare_best, collapse = '|')), moddirs)]purrr::map(moddirs_best, getcogVAST)Plot the centers of gravity by index
Small copepods in fall trending northeast similar to forage fish and aggregated survey species. Small copepods in spring trending north. Large copepods in fall trending west, similar to benthos. No significant change for Calfin, Euphausiids, Zooplankton volume.
purrr::map(moddirs_best, plotcog)## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
##
## [[12]]
##
## [[13]]
##
## [[14]]
##
## [[15]]
##
## [[16]]
##
## [[17]]
##
## [[18]]
Make the inputs for the SOE
from the “best” models for each season taxon combination.
Time series
SOEinputs <- function(d.name) { #, stratlook) {
# out of time to get this workign
# keep getting pmap error about x and y not the same source
# hard code it for SOE indices
#stratlook <- mget(stratlook)
# standard zoop strata for SOEs
stratlook <- data.frame(Stratum = c("Stratum_1",
"Stratum_2",
"Stratum_3",
"Stratum_4",
"Stratum_5",
"Stratum_6",
"Stratum_7"),
Region = c("AllEPU",
"her_sp",
"her_fa",
"MAB",
"GB",
"GOM",
"SS"))
modpath <- stringr::str_split(d.name, "/", simplify = TRUE)
modname <- modpath[length(modpath)]
outdir <- stringr::str_c(modpath[1:length(modpath)-1], collapse = "/")
infile <- file.path(d.name, "Index.csv")
season <- stringr::str_split(modname, "_")[[1]][2]
taxa <- stringr::str_split(modname, "_")[[1]][1]
outfile <- file.path(outdir, paste0(season, taxa, "index.rds"))
splitoutput <- read.csv(infile)
zoopindex <- splitoutput %>%
left_join(stratlook) %>%
dplyr::select(Time,
EPU = Region,
"Abundance Index Estimate" = Estimate,
"Abundance Index Estimate SE" = Std..Error.for.Estimate) %>%
tidyr::pivot_longer(c("Abundance Index Estimate", "Abundance Index Estimate SE"),
names_to = "Var", values_to = "Value") %>%
dplyr::filter(EPU %in% c("MAB", "GB", "GOM", "AllEPU")) %>%
dplyr::mutate(Units = "numbers per 100 cu m volume") %>%
dplyr::select(Time, Var, Value, EPU, Units)
zoopindex$Var <- stringr::str_c(stringr::str_to_title(season), stringr::str_to_title(taxa), zoopindex$Var, sep = " ")
saveRDS(zoopindex, outfile)
}
# # standard zoop strata for SOEs
# stratlook <- data.frame(Stratum = c("Stratum_1",
# "Stratum_2",
# "Stratum_3",
# "Stratum_4",
# "Stratum_5",
# "Stratum_6",
# "Stratum_7"),
# Region = c("AllEPU",
# "her_sp",
# "her_fa",
# "MAB",
# "GB",
# "GOM",
# "SS"))
#
# # larvarea zoop strata
# stratlook2 <- data.frame(Stratum = c("Stratum_1",
# "Stratum_2",
# "Stratum_3",
# "Stratum_4",
# "Stratum_5",
# "Stratum_6",
# "Stratum_7",
# "Stratum_8",
# "Stratum_9"),
# Region = c("AllEPU",
# "her_sp",
# "her_fa",
# "her_larv",
# "no_larv",
# "MAB",
# "GB",
# "GOM",
# "SS"))
#
#
# name <- dplyr::case_when(stringr::str_detect(moddirs_best, "larvarea") ~ "stratlook2",
# TRUE ~ "stratlook")
#
# list(d.name = moddirs_best,
# stratlook = name) |> purrr::pmap(SOEinputs)Center of gravity
SOEinputsCOG <- function(d.name) {
modpath <- stringr::str_split(d.name, "/", simplify = TRUE)
modname <- modpath[length(modpath)]
outdir <- stringr::str_c(modpath[1:length(modpath)-1], collapse = "/")
season <- stringr::str_split(modname, "_")[[1]][2]
taxa <- stringr::str_split(modname, "_")[[1]][1]
infile <- file.path(d.name, "cogout.rds")
cogout <- readRDS(infile)
outfile <- file.path(outdir, paste0(season, taxa, "cog.rds"))
groupcog <- as.data.frame(cogout$COG_Table) |>
dplyr::mutate(direction = ifelse(m==1, "Eastward", "Northward")) |>
dplyr::select("Time" = Year,
"Center of Gravity" = COG_hat,
"Center of Gravity SE" = SE,
direction) |>
tidyr::pivot_longer(c("Center of Gravity", "Center of Gravity SE"),
names_to = "Var", values_to = "Value") |>
#direction into Var
tidyr::unite(Var, direction:Var, sep = " ") |>
dplyr::mutate(Units = "km",
EPU = "ALLEPU") |>
dplyr::select(Time, Var, Value, EPU, Units)
groupcog$Var <- stringr::str_c(stringr::str_to_title(season), stringr::str_to_title(taxa), groupcog$Var, sep = " ")
#readr::write_csv(benthosindex, outfile)
saveRDS(groupcog, outfile)
}# we don't want the herring larval area run in the SOE, nor do we want herring larvae models
moddirs_SOE <- moddirs_best[!stringr::str_detect(moddirs_best, "larva")]
purrr::map(moddirs_SOE, SOEinputs)
purrr::map(moddirs_SOE, SOEinputsCOG)