From 444eced119a9cd415ce9bc1d91414338d716b91c Mon Sep 17 00:00:00 2001 From: kellijohnson-NOAA Date: Fri, 15 Apr 2022 05:55:16 -0700 Subject: [PATCH 01/10] fix(plot_qq): Uses `as.vector()` to remove class Some of the entries in `model[["data_list"]]` have the class "units" which makes performing operations on them difficult. Using `as.vector(...)` around the vectors with this class removes the errors that I was getting. See below ``` Error in Ops.units(model[["data_list"]][["b_i"]], 0) : both operands of the expression should be "units" objects ``` --- R/plot_qq.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/R/plot_qq.R b/R/plot_qq.R index 49803a2..b18000b 100644 --- a/R/plot_qq.R +++ b/R/plot_qq.R @@ -24,20 +24,20 @@ plot_qq <- function(outdir, model) { n_e <- model[["data_list"]][["n_c"]] e_i <- model[["data_list"]][["n_e"]] sigmaM <- (rep(1, n_e) %o% model[["Report"]][["SigmaM"]])[1, 1, 1] - Which <- which(model[["data_list"]][["b_i"]] > 0) + Which <- which(as.vector(model[["data_list"]][["b_i"]]) > 0) if (model[["data_list"]][["ObsModel_ez"]][1] == 1) { Qvals <- na.omit(mapply(plnorm, - q = model[["data_list"]][["b_i"]][Which], - meanlog = log(model[["data_list"]][["a_i"]][Which] * exp(model[["Report"]][["P2_iz"]][Which])) - + q = as.vector(model[["data_list"]][["b_i"]])[Which], + meanlog = log(as.vector(model[["data_list"]][["a_i"]])[Which] * exp(model[["Report"]][["P2_iz"]][Which])) - (sigmaM^2) / 2, MoreArgs = list(sdlog = sigmaM))) } if (model[["data_list"]][["ObsModel_ez"]][1] == 2) { Qvals <- na.omit(mapply(pgamma, - q = model[["data_list"]][["b_i"]][Which], + q = as.vector(model[["data_list"]][["b_i"]])[Which], scale = (sigmaM^2) * - (model[["data_list"]][["a_i"]][Which] * + (as.vector(model[["data_list"]][["a_i"]])[Which] * exp(model[["Report"]][["P2_iz"]][Which]) ), MoreArgs = list(shape = 1 / (sigmaM^2)))) From 263847ff24136af0bed8da75e7ca5c718e982ac6 Mon Sep 17 00:00:00 2001 From: kellijohnson-NOAA Date: Fri, 15 Apr 2022 06:11:31 -0700 Subject: [PATCH 02/10] feat: Saves model object if no standard errors are available Before, if there were no standard errors the model would continue instead of entering the early return. Errors/Warnings would arise because in the plotting of the index the standard errors are needed. We should not be doing diagnostics on models without inverted hessian matrices; so, now `VAST_do()` will exit early if the model did not converge rather than erroring out. --- R/VAST_do.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/VAST_do.R b/R/VAST_do.R index da47721..c3c3231 100644 --- a/R/VAST_do.R +++ b/R/VAST_do.R @@ -130,7 +130,8 @@ VAST_do <- function(Database, settings, conditiondir, compiledir, model_args = list(CompileDir = compiledir), silent = TRUE, run_model = TRUE), error = function(e) e) - if ("simpleError" %in% class(out)) { + if ("simpleError" %in% class(out) | + is.null(out$parameter_estimates$SD)) { save(out, file = file.path(conditiondir, "Error.RData")) save(list = ls(all.names = TRUE), file = file.path(conditiondir, "Save.RData")) return(out) From 6834506cde712992d93ad561b08ab900deb30fa2 Mon Sep 17 00:00:00 2001 From: kellijohnson-NOAA Date: Fri, 15 Apr 2022 06:13:50 -0700 Subject: [PATCH 03/10] fix: Provides `fit` to `FishStatsUtils::plot_biomass_index()` if `fit = NULL` which is the default, then the function doesn't work. This change provides the model output, which is labelled `out`, to the `fit` input argument. Closes #41 --- R/VAST_do.R | 1 + 1 file changed, 1 insertion(+) diff --git a/R/VAST_do.R b/R/VAST_do.R index c3c3231..8f08297 100644 --- a/R/VAST_do.R +++ b/R/VAST_do.R @@ -143,6 +143,7 @@ VAST_do <- function(Database, settings, conditiondir, compiledir, error = function(e) e) index <- suppressWarnings(FishStatsUtils::plot_biomass_index( + fit = out, DirName = file.path(conditiondir, .Platform$file.sep), TmbData = out$data_list, Sdreport = out$parameter_estimates$SD, use_biascorr = TRUE, From cb722e732a5c48606043a1b7f0bd94daf9c27825 Mon Sep 17 00:00:00 2001 From: kellijohnson-NOAA Date: Fri, 15 Apr 2022 06:16:57 -0700 Subject: [PATCH 04/10] fix(Table_for_SS3): Finds appropriate columns in indexdata Changes to VAST and FishStatsUtils that create the file Index rather than Table_for_SS3 led to column name changes. These legacy columns are easily recreated and saved similar to the old table in Table_for_SS3.csv to allow downstream code to work. This change was discussed in #39 with great changes made by @chantelwetzel-NOAA but we decided to go with the simpler change here rather than changing the column names of everything in VASTWestCoast which would also lead to code changes that would be needed for assessment purposes. --- R/VAST_do.R | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/R/VAST_do.R b/R/VAST_do.R index 8f08297..c75e6b7 100644 --- a/R/VAST_do.R +++ b/R/VAST_do.R @@ -156,8 +156,14 @@ VAST_do <- function(Database, settings, conditiondir, compiledir, parameter_estimates = out$parameter_estimates, savedir = conditiondir) fileindex <- file.path(conditiondir, "Table_for_SS3.csv") - indexdata <- utils::read.csv(fileindex) - indexdata[,"Year"] <- out$year_labels[indexdata[,"Year"]] + indexdata <- data.frame( + Year = out$year_labels[index[["Table"]][, "Time"]], + Unit = index[["Table"]][, "Units"], + Fleet = index[["Table"]][, "Stratum"], + Estimate_metric_tons = index[["Table"]][, "Estimate"], + SD_log = index[["Table"]][, "Std. Error for Estimate"], + SD_mt = index[["Table"]][, "Std. Error for ln(Estimate)"] + ) utils::write.csv(x = indexdata, file = fileindex, row.names = FALSE) plot_ss(file.in = fileindex, lab.survey = survey, From d055d6dbb3aafca47505b19c541ddf3893fcf959 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Fri, 15 Apr 2022 07:48:51 -0700 Subject: [PATCH 05/10] refactor: change plot area colors --- R/plot_ss.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/plot_ss.R b/R/plot_ss.R index b940f5f..6a8c8d0 100644 --- a/R/plot_ss.R +++ b/R/plot_ss.R @@ -49,19 +49,19 @@ plot_ss <- function(file.in = "Table_for_SS3.csv", savefile, savefile <- normalizePath(savefile, mustWork = FALSE) dir.create(dirname(savefile), recursive = TRUE, showWarnings = FALSE) - colors <- c("black", "blue", "darkorchid1", "red") + colors <- c("black", "red", "darkorchid1", "blue") pch.vec <- c(21, 22, 23, 24) cex.vec <- c(1.6, 1.4, 1.4, 1.4) cex.lab <- 2.1 cex.axis <- 1.6 - colors <- c("black", "blue", "darkorchid1", "red") + if (length(strat) > 4) { stop("plot_ss is not designed to work with more than 4 strata.", "\nPlease contact the maintainer(s) to request additional", "\nfunctionality here or limit the number of strata in the csv.") } if ("north" %in% tolower(strat) && length(strat) == 3) { - colors <- c("black", "blue", "red") + colors <- c("black", "red", "blue") } grDevices::png(filename = savefile, From c651793415ebb052f2deabf7ff7f9bfaf9958fcb Mon Sep 17 00:00:00 2001 From: kellijohnson-NOAA Date: Fri, 15 Apr 2022 08:57:43 -0700 Subject: [PATCH 06/10] doc(zone): Adds comment on why we do not use zone --- R/VAST_do.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/VAST_do.R b/R/VAST_do.R index c75e6b7..5b22722 100644 --- a/R/VAST_do.R +++ b/R/VAST_do.R @@ -74,6 +74,8 @@ VAST_do <- function(Database, settings, conditiondir, compiledir, fine_scale = settings[["fine_scale"]], strata.limits = settings[["strata"]], #zone = NA, #default + # automatically detecting zone in `FishStatsUtils::make_extrapolation_info()` + # DO NOT allow users to control just use default FieldConfig = settings[["FieldConfig"]], RhoConfig = settings[["RhoConfig"]], OverdispersionConfig = settings[["overdispersion"]], From c359cdeca5ec82e809b1f559b443e451d251069b Mon Sep 17 00:00:00 2001 From: kellijohnson-NOAA Date: Fri, 15 Apr 2022 13:01:35 -0700 Subject: [PATCH 07/10] feat(use_anisotropy): Allows users to change TRUE use_anisotropy = TRUE is the default in make_settings(), which was previously used by VASTWestCoast, now users can put FALSE in their get_settings() input list if they want to chagne it. This turns off the estimation of the H matrix that allows for spatial rotation, i.e., rotation of the axis to be slanted rather than perpindicular to the coast. --- R/VAST_do.R | 2 +- R/get_settings.R | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/R/VAST_do.R b/R/VAST_do.R index 5b22722..0261a2b 100644 --- a/R/VAST_do.R +++ b/R/VAST_do.R @@ -83,7 +83,7 @@ VAST_do <- function(Database, settings, conditiondir, compiledir, bias.correct = TRUE, #calculate derived quantities of interest, no harm #Options = ,#default is more than original VASTWestCoast - #use_anisotropy = TRUE, #default + use_anisotropy = settings[["use_anisotropy"]], Version = settings[["version"]], #treat_nonencounter_as_zero = TRUE, #default #VamConfig = c(Method = 0, Rank = 0, Timing = 0), #default diff --git a/R/get_settings.R b/R/get_settings.R index e7e1203..a205fc8 100644 --- a/R/get_settings.R +++ b/R/get_settings.R @@ -32,7 +32,9 @@ get_settings <- function(settings = NULL, verbose = FALSE) { "field" = NULL, "rho" = NULL, "fine_scale" = TRUE, - "overdispersion" = NULL) + "overdispersion" = NULL, + "use_anisotropy" = TRUE + ) need <- !names(Settings_all) %in% names(settings) if (verbose) { message("Adding the following objects to settings:\n", From e4eeae9b50029094936da41b5a519baaa7175d34 Mon Sep 17 00:00:00 2001 From: kellijohnson-NOAA Date: Fri, 15 Apr 2022 14:15:15 -0700 Subject: [PATCH 08/10] Adds Tweedie closes #43 --- R/VAST_spp.R | 5 +++-- R/get_settings.R | 5 +++++ R/plot_qq.R | 4 ++++ R/summary_nwfsc.R | 4 +++- man/VAST_spp.Rd | 2 +- 5 files changed, 16 insertions(+), 4 deletions(-) diff --git a/R/VAST_spp.R b/R/VAST_spp.R index 21a6c45..039fd00 100644 --- a/R/VAST_spp.R +++ b/R/VAST_spp.R @@ -85,7 +85,7 @@ #' } #' VAST_spp <- function(dir, species, - dist = c("lognormal", "gamma")) { + dist = c("lognormal", "gamma", "tweedie")) { #### Get species, survey, and strata info info <- nwfscSurvey::GetSpp.fn(species) @@ -107,7 +107,8 @@ VAST_spp <- function(dir, species, for (obs in dist) { obs_model <- switch(obs, lognormal = c(1, 0), - gamma = c(2, 0) + gamma = c(2, 0), + tweedie = c(10, 2) ) modeldir <- file.path(sppdir, paste(survey, obs, sep = "_")) settings <- list( diff --git a/R/get_settings.R b/R/get_settings.R index a205fc8..6e67c5a 100644 --- a/R/get_settings.R +++ b/R/get_settings.R @@ -69,6 +69,11 @@ get_settings <- function(settings = NULL, verbose = FALSE) { default = c(Beta1 = 0, Beta2 = 0, Epsilon1 = 0, Epsilon2 = 0)) Settings_all$FieldConfig <- get_settings_single(Settings_all$field, default = c(Omega1 = 1, Epsilon1 = 1, Omega2 = 1, Epsilon2 = 1)) + if (Settings_all[["ObsModelcondition"]][1] == 10) { + # Following items are specific to the Tweedie distribution #43 + Settings_all[["FieldConfig"]][c("Omega1", "Epsilon1")] <- 0 + Settings_all[["RhoConfig"]][c("Beta1")] <- 3 + } # Overdispersion if (is.null(Settings_all[["overdispersion"]])) { diff --git a/R/plot_qq.R b/R/plot_qq.R index b18000b..33c922f 100644 --- a/R/plot_qq.R +++ b/R/plot_qq.R @@ -42,6 +42,10 @@ plot_qq <- function(outdir, model) { ), MoreArgs = list(shape = 1 / (sigmaM^2)))) } + if (model[["data_list"]][["ObsModel_ez"]][1] == 10) { + message("QQ plot is not yet implemented for the Tweedie distribution. See #45.") + return(invisible()) + } grDevices::png(file.path(outdir, "VASTWestCoast_QQ.png"), height = 18, width = 18, units = "cm", res = 300) diff --git a/R/summary_nwfsc.R b/R/summary_nwfsc.R index d2f8937..9f81cb4 100644 --- a/R/summary_nwfsc.R +++ b/R/summary_nwfsc.R @@ -42,7 +42,9 @@ summary_nwfsc <- function(obj, parameter_estimates, savedir = NULL) { "Poisson-link", switch(as.character(obj$env$data$ObsModel[1]), "1" = "Lognormal", - "2" = "Gamma"))) + "2" = "Gamma", + "10" = "Tweedie")) + ) TableA[6, ] <- c("Spatial effect for encounter probability", ifelse(obj$env$data$FieldConfig[1, 2] < 0, "No", "Yes")) TableA[7, ] <- c("Spatio-temporal effect for encounter probability", diff --git a/man/VAST_spp.Rd b/man/VAST_spp.Rd index 8e187a5..cce5590 100644 --- a/man/VAST_spp.Rd +++ b/man/VAST_spp.Rd @@ -4,7 +4,7 @@ \alias{VAST_spp} \title{Run \pkg{VAST} models for a species sampled by west coast surveys} \usage{ -VAST_spp(dir, species, dist = c("lognormal", "gamma")) +VAST_spp(dir, species, dist = c("lognormal", "gamma", "tweedie")) } \arguments{ \item{dir}{An existing directory where you want the new folder that is created From 44530d0ff82248923b31b5ba88ad0966902283cb Mon Sep 17 00:00:00 2001 From: kellijohnson-NOAA Date: Thu, 21 Apr 2022 13:50:56 -0700 Subject: [PATCH 09/10] fix(Table_for_SS3.csv): Use mt Close #47 thanks to @chantelwetzel-noaa Close #48 --- R/VAST_do.R | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/R/VAST_do.R b/R/VAST_do.R index 0261a2b..5b502bd 100644 --- a/R/VAST_do.R +++ b/R/VAST_do.R @@ -162,9 +162,10 @@ VAST_do <- function(Database, settings, conditiondir, compiledir, Year = out$year_labels[index[["Table"]][, "Time"]], Unit = index[["Table"]][, "Units"], Fleet = index[["Table"]][, "Stratum"], - Estimate_metric_tons = index[["Table"]][, "Estimate"], - SD_log = index[["Table"]][, "Std. Error for Estimate"], - SD_mt = index[["Table"]][, "Std. Error for ln(Estimate)"] + # Go from kg to mt to keep backwards compatibility with VASTWestCoast + Estimate_metric_tons = index[["Table"]][, "Estimate"] / 1000, + SD_log = index[["Table"]][, "Std. Error for ln(Estimate)"], + SD_mt = index[["Table"]][, "Std. Error for Estimate"] / 1000 ) utils::write.csv(x = indexdata, file = fileindex, row.names = FALSE) plot_ss(file.in = fileindex, From 79fcb086004a97d517ef989d0246a4f341e4f84a Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Fri, 13 May 2022 06:46:14 -0700 Subject: [PATCH 10/10] fix(Table_for_SS4.csv): change output units to mt from kg --- R/VAST_do.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/VAST_do.R b/R/VAST_do.R index 5b502bd..5b60167 100644 --- a/R/VAST_do.R +++ b/R/VAST_do.R @@ -160,7 +160,8 @@ VAST_do <- function(Database, settings, conditiondir, compiledir, fileindex <- file.path(conditiondir, "Table_for_SS3.csv") indexdata <- data.frame( Year = out$year_labels[index[["Table"]][, "Time"]], - Unit = index[["Table"]][, "Units"], + # Change units to mt since converting from kg to mt below + Unit = "mt", Fleet = index[["Table"]][, "Stratum"], # Go from kg to mt to keep backwards compatibility with VASTWestCoast Estimate_metric_tons = index[["Table"]][, "Estimate"] / 1000,