diff --git a/.Rbuildignore b/.Rbuildignore index 2c126f5f9..44d896c68 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -2,6 +2,7 @@ ^\.Rproj\.user$ ^\.github$ ^\.history$ +^\.agents$ CMakeLists.txt tests/gtest tests/test_plan diff --git a/DESCRIPTION b/DESCRIPTION index ee12e15f6..6ea9c96a1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: FIMS Title: The Fisheries Integrated Modeling System -Version: 0.9.4.9000 +Version: 0.10.0 Authors@R: c( person(c("Kelli", "F."), "Johnson", , "kelli.johnson@noaa.gov", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-5149-451X")), @@ -72,7 +72,7 @@ URL: https://github.com/noaa-fims/fims, https://noaa-fims.github.io, https://noaa-fims.github.io/FIMS BugReports: https://github.com/noaa-fims/fims/issues Depends: - R (>= 4.1.0) + R (>= 4.2.0) Imports: cli, dplyr, diff --git a/NAMESPACE b/NAMESPACE index bf937c7f8..82199505b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -60,9 +60,9 @@ export(get_parameter_names) export(get_random) export(get_random_names) export(get_report) +export(get_run_time) export(get_sdreport) export(get_start_year) -export(get_timing) export(get_version) export(glance) export(initialize_data_distribution) diff --git a/NEWS.md b/NEWS.md index f3955ea5e..968e0068a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,7 @@ # FIMS 0.10.0 -* +* Use setup_default_parameters() to set up a model +* Specify uncertainty using a formula with a distribution and arguments # FIMS 0.9.4 diff --git a/R/check_convergence.R b/R/check_convergence.R index f6b7c7c90..df99824b4 100644 --- a/R/check_convergence.R +++ b/R/check_convergence.R @@ -90,7 +90,7 @@ check_mle_convergence <- function(input, obj, opt, maxgrad) { obj = obj, opt = opt, sdreport = list(), - timing = c( + run_time = c( time_optimization = as.difftime(0, units = "secs"), time_sdreport = as.difftime(0, units = "secs"), time_total = as.difftime(0, units = "secs") @@ -151,7 +151,7 @@ check_sdreport_convergence <- function(input, obj, opt, sdreport) { obj = obj, opt = opt, sdreport = sdreport, - timing = c( + run_time = c( time_optimization = as.difftime(0, units = "secs"), time_sdreport = as.difftime(0, units = "secs"), time_total = as.difftime(0, units = "secs") diff --git a/R/data_.R b/R/data_.R index b5d77a05a..f4a53728c 100644 --- a/R/data_.R +++ b/R/data_.R @@ -31,12 +31,12 @@ #' `landings`, and `weight_at_age` data. `number` or `proportion` are each #' viable units for the composition data, where the former is the preferred #' unit of measurement.} -#' \item{uncertainty}{A real value providing a measurement of uncertainty -#' for value. For landings and indices of abundance this should be -#' the standard deviation of the logged observations if you are using the -#' lognormal distribution to fit your data. For composition data it will -#' be your input sample size. +#' \item{uncertainty}{A right-handed formula specifying the distributional +#' assumptions for the entry in `value`. See [FIMSFrame()] for more +#' information about this column. #' } #' } #' @source \url{www.github.com/NOAA-FIMS/Age_Structured_Stock_Assessment_Model_Comparison} +#' @seealso +#' * [FIMSFrame()] "data_big" diff --git a/R/distribution_formulas.R b/R/distribution_formulas.R index a29fbcb09..e2ab03aa6 100644 --- a/R/distribution_formulas.R +++ b/R/distribution_formulas.R @@ -25,26 +25,10 @@ check_distribution_validity <- function(args) { # Separate objects from args family <- args[["family"]] sd <- args[["sd"]] - # Optional argument data_type - data_type <- args[["data_type"]] check_present <- purrr::map_vec(list("family" = family, "sd" = sd), is.null) - # Set up global rules - # FIXME: Move this to a data item in the package so it can be used everywhere - # Could do a call to all data objects in the package and get unique types that - # are available - data_type_names <- c("landings", "index", "agecomp", "lengthcomp") - if (is.null(data_type)) { - available_distributions <- c("lognormal", "gaussian") - } else { - available_distributions <- switch( - EXPR = ifelse(grepl("comp", data_type), "composition", data_type), - "landings" = c("lognormal", "gaussian"), - "index" = c("lognormal", "gaussian"), - "composition" = c("multinomial"), - "unavailable data type" - ) - } + # Only process distributions are currently validated here. + available_distributions <- c("lognormal", "gaussian") elements_of_sd <- c("value", "estimation_type") # Start a bulleted list of errors and add to it in each if statement @@ -71,38 +55,18 @@ check_distribution_validity <- function(args) { e.g., `family = gaussian()`, instead of {class(family)}." ) } else { - if ( - !(family[["family"]] %in% available_distributions) || - "unavailable data type" %in% available_distributions - ) { - ifelse_type <- ifelse( - is.null(data_type), - "distribution", - paste(data_type, "data") - ) + family_name <- family[["family"]] + if (!(family_name %in% available_distributions)) { abort_bullets <- c( abort_bullets, "x" = "FIMS currently does not allow the family to be - {.code {family[['family']]}}.", - "i" = "The families available for this {ifelse_type} are + {.code {family_name}}.", + "i" = "The families available for process distributions are {.code {available_distributions}}." ) } } - # Checks related to the type of data - if (!is.null(data_type)) { - if (!(data_type %in% data_type_names)) { - abort_bullets <- c( - abort_bullets, - "x" = "The specified {.var data_type} of {.var {data_type}} is not - available.", - "i" = "Allowed values for {.var data_type} are - {.code {data_type_names}}." - ) - } - } - # Checks related to standard deviation # Check if sd has both elements and if yes, then go onto the else statement # for major checks @@ -139,25 +103,6 @@ check_distribution_validity <- function(args) { } } - # Check dimensions for data distributions: sd must be either length 1 (scalar) - # or match data length when data_type is landings or index - if (!is.null(data_type) && !is.null(args[["module"]])) { - module <- args[["module"]] - if (data_type == "landings" | data_type == "index") { - n_obs <- module$n_years$get() - - if (length(sd[["value"]]) > 1 && length(sd[["value"]]) != n_obs) { - abort_bullets <- c( - abort_bullets, - "x" = "The size of {.var log_sd} does not match the size of observed data for {data_type}.", - "i" = "The {.var log_sd} vector is of size {length(sd[['value']])}.", - "i" = "The observed {data_type} data vector is of size {n_obs}.", - "i" = "Either provide a single {.var log_sd} value (scalar) or a vector matching the data length." - ) - } - } - } - # Return error messages if more than just the default is present if (length(abort_bullets) == 1) { invisible(TRUE) @@ -166,52 +111,6 @@ check_distribution_validity <- function(args) { } } -#' Return name of expected value -#' -#' The combination of data type, family, and link lead to a specific name for -#' the expected value within the code base. This function looks at the -#' combination of these three objects and specifies the appropriate string for -#' its name going forward. -#' @inheritParams initialize_data_distribution -#' @noRd -#' @return -#' A string specifying the name of the expected value. -#' -get_expected_name <- function(family, data_type) { - # TODO: Think about if the name of the expected value should change based on - # the link or if it should stay the same? Keeping track of different names in - # the code base might be too complex for the output as well - family_string <- family[["family"]] - link_string <- family[["link"]] - expected_name <- dplyr::case_when( - data_type == "landings" && - grepl("lognormal|gaussian", family_string) && - link_string == "log" ~ "log_landings_expected", - data_type == "landings" && - grepl("lognormal|gaussian", family_string) && - link_string == "identity" ~ "landings_expected", - data_type == "index" && - grepl("lognormal|gaussian", family_string) && - link_string == "log" ~ "log_index_expected", - data_type == "index" && - grepl("lognormal|gaussian", family_string) && - link_string == "identity" ~ "index_expected", - grepl("agecomp", data_type) ~ "agecomp_proportion", - grepl("lengthcomp", data_type) ~ "lengthcomp_proportion", - ) - # Check combination of entries was okay and led to valid name - if (is.na(expected_name)) { - cli::cli_abort(c( - "x" = "The combination of data type, family, and link are incompatible in - some way.", - "i" = "{.var data_type} is {.var {data_type}}.", - "i" = "The family is {.var {family_string}}.", - "i" = "The link is {.var {link_string}}." - )) - } - return(expected_name) -} - #' Set up a new distribution for a data type or a process #' #' Use [methods::new()] to set up a distribution within an existing module with @@ -223,26 +122,14 @@ get_expected_name <- function(family, data_type) { #' `initialize_process_distribution()`. #' @param module An identifier to a C++ fleet module that is linked to the data #' of interest. -#' @param family A description of the error distribution and link function to -#' be used in the model. The argument takes a family class, e.g., -#' `stats::gaussian(link = "identity")`. -#' @param sd A list of length two. The first entry is named `"value"` and it -#' stores the initial values (scalar or vector) for the relevant standard -#' deviations. The default is `value = 1`. The second entry is named -#' `"estimation_type"` and it stores a vector of booleans (default = -#' "constant") is a string indicating whether or not standard deviation is -#' estimated as a fixed effect or held constant. If `"value"` is a vector and -#' `"estimation_type"` is a scalar, the single value specified -#' `"estimation_type"` value will be repeated to match the length of `value`. -#' Otherwise, the dimensions of the two must match. #' @param data_type A string specifying the type of data that the #' distribution will be fit to. Allowable types include #' `r glue::glue_collapse(sprintf('"%s"', eval(formals(initialize_data_distribution)[["data_type"]])), sep = ", ", last = ", and ")` #' and the default is #' `r eval(formals(initialize_data_distribution)[["data_type"]])[1]`. -#' @param par A string specifying the parameter name the distribution applies -#' to. Parameters must be members of the specified module. Use -#' `methods::show(module)` to obtain names of parameters within the module. +#' @param uncertainty A vector of strings specifying formulas for each data +#' point. See [FIMSFrame()] for more information on what the formula should +#' look like. #' @return #' A reference class. is returned. Use [methods::show()] to view the various #' Rcpp class fields, methods, and documentation. @@ -255,12 +142,11 @@ get_expected_name <- function(family, data_type) { #' # Create a new fleet module #' fleet <- methods::new(Fleet) #' # Create a distribution for the fleet module +#' sd_log <- rep(sqrt(log(0.01^2 + 1)), n_years) #' fleet_distribution <- initialize_data_distribution( #' module = fishing_fleet, -#' family = lognormal(link = "log"), -#' sd = list( -#' value = rep(sqrt(log(0.01^2 + 1)), n_years), -#' estimation_type = rep("constant", n_years) # Can be a single "constant" +#' uncertainty = glue::glue( +#' "~dlnorm(meanlog = log_index_expected, sdlog = {sd_log})" #' ), #' data_type = "index" #' ) @@ -280,78 +166,47 @@ get_expected_name <- function(family, data_type) { #' } initialize_data_distribution <- function( module, - family = NULL, - # Create a tibble with value and estimation_type column for sd - sd = tibble::tibble( - value = 1, - estimation_type = "constant" - ), - # FIXME: Move this argument to second to match where par is in - # initialize_process_distribution - data_type = c("landings", "index", "agecomp", "lengthcomp") + data_type = c("landings", "index", "age_comp", "length_comp"), + uncertainty ) { data_type <- rlang::arg_match(data_type) - # FIXME: Make the available families a data object - # Could also make the matrix of distributions available per type as a - # data frame where the check could use the stored object. - - # validity check on user input - args <- list( - family = family, - sd = sd, - data_type = data_type, - module = module - ) - check_distribution_validity(args) - - # assign name of observed data based on data_type - obs_id_name <- glue::glue("observed_{data_type}_data_id") + uncertainty_split <- purrr::map( + .x = as.list(uncertainty), + .f = parse_data_distribution + ) |> + transpose_data_distribution() # Set up distribution based on `family` argument` - if (family[["family"]] == "lognormal") { - # create new Rcpp module - new_module <- methods::new(DlnormDistribution) - - # populate logged standard deviation parameter with log of input - # Using resize() and then assigning value to each element of log_sd directly - # is correct, as creating a new VariableVector for log_sd here would - # trigger an error in integration tests with wrappers. - new_module$log_sd[] <- log(sd[["value"]]) - new_module$log_sd$set_estimation_types(sd[["estimation_type"]]) - } - - if (family[["family"]] == "gaussian") { - # create new Rcpp module - new_module <- methods::new(DnormDistribution) - - # populate logged standard deviation parameter with log of input - new_module$log_sd[] <- log(sd[["value"]]) - new_module$log_sd$set_estimation_types(sd[["estimation_type"]]) - } - - if (family[["family"]] == "multinomial") { - # create new Rcpp module - new_module <- methods::new(DmultinomDistribution) + new_module <- methods::new(get( + gsub("^d(.+)$", "D\\1Distribution", uncertainty_split[["family"]][[1]]) + )) + if ("log_sd" %in% names(new_module)) { + new_module$log_sd[] <- log( + unlist(uncertainty_split[[grep("sd", names(uncertainty_split))]]) + ) + new_module$log_sd$set_estimation_types("constant") } # setup link to observed data - if (data_type == "landings") { - new_module$set_observed_data(module$GetObservedLandingsDataID()) - } - if (data_type == "index") { - new_module$set_observed_data(module$GetObservedIndexDataID()) - } - if (data_type == "agecomp") { - new_module$set_observed_data(module$GetObservedAgeCompDataID()) - } - if (data_type == "lengthcomp") { - new_module$set_observed_data(module$GetObservedLengthCompDataID()) + data_id_name <- glue::glue("GetObserved{snake_to_pascal(data_type)}DataID") + getter_function <- tryCatch( + do.call("$", list(module, data_id_name)), + error = function(e) NULL + ) + if (is.null(getter_function)) { + cli::cli_abort(c( + x = "Could not find observed data name, {.code {data_id_name}} in the + module", + i = "Check {.var data_type} and the module class.", + i = "Options in the module are {grep('Get', names(module), value = TRUE)}" + )) } - - # set name of expected values - expected <- get_expected_name(family, data_type) + new_module$set_observed_data(getter_function()) # setup link to expected values - new_module$set_distribution_links("data", module$field(expected)$get_id()) + new_module$set_distribution_links( + "data", + module$field(uncertainty_split[["link"]][[1]])$get_id() + ) return(new_module) } diff --git a/R/fims_distributions.R b/R/fims_distributions.R new file mode 100644 index 000000000..c6e3866e9 --- /dev/null +++ b/R/fims_distributions.R @@ -0,0 +1,11 @@ +#' Distribution parameter lookup table +#' +#' A lookup table that maps each supported distribution family to the field +#' name used to link observed data, the random-value generator, and the +#' remaining parameter names for that family. +#' +#' @format A data frame with columns `family`, `sample_function`, `central_parameter`, and +#' `other_parameters`. `sample_function` is a character column with function +#' names. +#' @keywords datasets +"fims_distributions" diff --git a/R/fimsfit.R b/R/fimsfit.R index 388299269..48ceb3adb 100644 --- a/R/fimsfit.R +++ b/R/fimsfit.R @@ -31,7 +31,7 @@ methods::setClass( report = "list", sdreport = "sdreportOrList", number_of_parameters = "integer", - timing = "difftime", + run_time = "difftime", version = "package_version", model_output = "character" ) @@ -59,7 +59,7 @@ methods::setMethod( f = "print", signature = "FIMSFit", definition = function(x) { - rt <- as.numeric(x@timing[["time_total"]], units = "secs") + rt <- as.numeric(x@run_time[["time_total"]], units = "secs") ru <- "seconds" if (rt > 60 * 60 * 24) { rt <- rt / (60 * 60 * 24) @@ -285,15 +285,15 @@ methods::setMethod( ) #' @return -#' [get_timing()] returns the amount of time it took to run the model in +#' [get_run_time()] returns the amount of time it took to run the model in #' seconds as a `difftime` object. #' @export #' @rdname get_FIMSFit #' @keywords fit_fims -methods::setGeneric("get_timing", function(x) standardGeneric("get_timing")) +methods::setGeneric("get_run_time", function(x) standardGeneric("get_run_time")) #' @rdname get_FIMSFit #' @keywords fit_fims -methods::setMethod("get_timing", "FIMSFit", function(x) x@timing) +methods::setMethod("get_run_time", "FIMSFit", function(x) x@run_time) #' @return #' [get_version()] returns the `package_version` of FIMS that was used to fit @@ -370,7 +370,7 @@ is.FIMSFit <- function(x) { #' [stats::nlminb()], used to fit a TMB model. #' @param sdreport An object of the `sdreport` class as returned from #' [TMB::sdreport()]. -#' @param timing A vector of at least length one, where all entries are of the +#' @param run_time A vector of at least length one, where all entries are of the #' `timediff` class and at least one is named "time_total". This information #' is available in [fit_fims()] and added to this argument internally but if #' you are a power user you can calculate the time it took to run your model @@ -403,7 +403,7 @@ is.FIMSFit <- function(x) { #' An object with the `sdreport` class containing the output from #' `TMB::sdreport(obj)`. #' } -#' \item{\code{timing}:}{ +#' \item{\code{run_time}:}{ #' The length of time it took to run the model if it was optimized. #' } #' \item{\code{version}:}{ @@ -422,7 +422,7 @@ FIMSFit <- function( obj, opt = list(), sdreport = list(), - timing = c("time_total" = as.difftime(0, units = "secs")), + run_time = c("time_total" = as.difftime(0, units = "secs")), version = utils::packageVersion("FIMS") ) { # Determine the number of parameters @@ -480,7 +480,7 @@ FIMSFit <- function( report = report, sdreport = sdreport, number_of_parameters = number_of_parameters, - timing = timing, + run_time = run_time, version = version, model_output = model_output ) @@ -562,7 +562,7 @@ fit_fims <- function(input, initial_fit <- FIMSFit( input = input, obj = obj, - timing = c("time_total" = as.difftime(0, units = "secs")) + run_time = c("time_total" = as.difftime(0, units = "secs")) ) return(initial_fit) } @@ -585,7 +585,7 @@ fit_fims <- function(input, obj = obj, opt = failed_nlminb_object[["opt"]], sdreport = list(), - timing = failed_nlminb_object[["timing"]] + run_time = failed_nlminb_object[["run_time"]] ) return(failed_fit) } @@ -619,7 +619,7 @@ fit_fims <- function(input, obj = obj, opt = failed_nlminb_object[["opt"]], sdreport = list(), - timing = failed_nlminb_object[["timing"]] + run_time = failed_nlminb_object[["run_time"]] ) return(failed_fit) } @@ -651,7 +651,7 @@ fit_fims <- function(input, time_sdreport <- as.difftime(0, units = "secs") } - timing <- c( + run_time <- c( time_optimization = time_optimization, time_sdreport = time_sdreport, time_total = Sys.time() - t0 @@ -661,7 +661,7 @@ fit_fims <- function(input, obj = obj, opt = opt, sdreport = sdreport, - timing = timing + run_time = run_time ) print(fit) if (!is.null(filename)) { @@ -739,7 +739,7 @@ return_failed_nlminb <- function(object) { # Construct a fallback optimizer result with consistent structure, i.e., # convergence = 1L indicates non-convergence return(list( - timing = c( + run_time = c( time_optimization = as.difftime(0, units = "secs"), time_sdreport = as.difftime(0, units = "secs"), time_total = as.difftime(0, units = "secs") diff --git a/R/fimsframe.R b/R/fimsframe.R index 422460a7e..b6f375e8a 100644 --- a/R/fimsframe.R +++ b/R/fimsframe.R @@ -399,7 +399,23 @@ methods::setMethod( .data[["type"]] == "age_comp", .data[["fleet"]] %in% .env$fleet ) |> - dplyr::pull(.data[["value"]]) + dplyr::mutate( + sample_size = ifelse( + .data$value == -999, + 1, + purrr::map_dbl( + purrr::map(.f = parse_data_distribution, .x = .data$uncertainty), + list("size"), + .default = NA_real_ + ) + ), + probabilities = ifelse( + .data$unit == "number", + .data$value, + .data$value * .data$sample_size + ) + ) |> + dplyr::pull(.data[["probabilities"]]) } ) #' @rdname model_ @@ -438,7 +454,23 @@ methods::setMethod( .data[["type"]] == "length_comp", .data[["fleet"]] %in% .env$fleet ) |> - dplyr::pull(.data[["value"]]) + dplyr::mutate( + sample_size = ifelse( + .data$value == -999, + 1, + purrr::map_dbl( + purrr::map(.f = parse_data_distribution, .x = .data$uncertainty), + list("size"), + .default = NA_real_ + ) + ), + probabilities = ifelse( + .data$unit == "number", + .data$value, + .data$value * .data$sample_size + ) + ) |> + dplyr::pull(.data[["probabilities"]]) } ) #' @rdname model_ @@ -858,7 +890,8 @@ validate_dimension_of_conversion <- function(data, n_groups, n_timings) { #' @details #' ## data #' The input data are both sorted (see the section below on sorting) and -#' expanded before returning them in the data slot. +#' expanded to include -999 values for all missing rows before returning them +#' in the data slot. #' ### Ages #' Currently, ages must be integers, i.e., FIMS cannot accommodate numeric ages #' like age 1.5 but we hope that this is something that we will be able to @@ -870,6 +903,28 @@ validate_dimension_of_conversion <- function(data, n_groups, n_timings) { #' age in the model. For example, you cannot bin fish into bins that span #' multiple years, you must have age-2, age-3, and age-4 not just age-2 and #' age-4 fish in your composition data. +#' ### Uncertainty +#' Uncertainty information for your data contains information for fitting the +#' model and for creating bootstrapped data sets. Right-handed formulas are used +#' to specify the distribution your data is assumed to follow, the derived +#' quantity in the model that the data is assumed to represent, and parameters +#' associated with the assumed distribution. For example, landings data can be +#' distributed using a normal or lognormal distribution and if you use a +#' normal distribution the value is assumed to represent `landings_expected` but +#' if you use a lognormal distribution the value is assumed to represent +#' `log_landings_expected`. Therefore, for the former you would specify +#' `"~dnorm(mean = landings_expected, sd = 0.01)"` and for the later you would +#' specify `"~dlnorm(meanlog = log_landings_expected, sdlog = 0.01)"`. Where, +#' the input value for sd or sdlog is the standard deviation of the mean value +#' or standard deviation of the log of the mean value, respectively. The +#' distribution names and argument names match the names used in R. Use +#' [args()], e.g., `args(dnorm)`, to get information on the argument names that +#' are required for your distribution. And, use `?distributions` to get a list +#' of distributions that have density functions in base R. To get information on +#' available derived quantities to set the central tendency of the distribution +#' to, you can create a fleet object and inspect the available names, e.g., +#' `names(methods::new(Fleet))`, or you can navigate to the Public Attributes +#' section of the doxygen documentation for Fleet by running `?Fleet`. #' ### Sorting #' It is important that the order of the rows in the data are correct but it is #' not expected that the user will do this. Instead, the returned data are diff --git a/R/initialize_modules.R b/R/initialize_modules.R index 7f2718274..3f3d75d2b 100644 --- a/R/initialize_modules.R +++ b/R/initialize_modules.R @@ -518,18 +518,7 @@ initialize_comp <- function(data, )) } - model_data <- comp_data * - get_data(data) |> - dplyr::filter( - .data$fleet == .env$fleet, - .data$type == comp[["name"]] - ) |> - dplyr::mutate( - valid_n = ifelse(.data$value == -999, 1, .data$uncertainty) - ) |> - dplyr::pull(.data$valid_n) - - if (length(model_data) != get_n_years(data) * get_function(data)) { + if (length(comp_data) != get_n_years(data) * get_function(data)) { bad_data_years <- get_data(data) |> dplyr::filter( .data$fleet == .env$fleet, @@ -543,12 +532,12 @@ initialize_comp <- function(data, "The length of the `{comp[['name']]}`-composition data for fleet `{fleet}` does not match the expected dimensions.", i = "Expected length: {get_n_years(data) * get_function(data)}", - i = "Actual length: {length(model_data)}", - i = "Number of -999 values: {sum(model_data == -999)}", + i = "Actual length: {length(comp_data)}", + i = "Number of -999 values: {sum(comp_data == -999)}", i = "Dates with invalid data: {bad_data_years}" )) } - module[[comp[["comp_data_field"]]]][] <- model_data + module[[comp[["comp_data_field"]]]][] <- comp_data return(module) } @@ -739,61 +728,55 @@ initialize_fims <- function(parameters, data) { if ("index" %in% fleet_types) { fleet_index_distribution[[i]] <- initialize_data_distribution( + data_type = "index", module = fleet[[i]], # TODO: need to update family and match options from the distribution # column from the parameters tibble - family = lognormal(link = "log"), - sd = get_data(data) |> + uncertainty = get_data(data) |> dplyr::filter( .data$fleet == .env$fleets[i] & .data$type == "index" ) |> - dplyr::select(-value) |> - dplyr::mutate( - estimation_type = "constant", - value = .data$uncertainty - ), - data_type = "index" + dplyr::pull(dplyr::all_of("uncertainty")) ) } if ("landings" %in% fleet_types) { fleet_landings_distribution[[i]] <- initialize_data_distribution( module = fleet[[i]], - # TODO: need to update family and match options from the distribution - # column from the parameters tibble - family = lognormal(link = "log"), - sd = get_data(data) |> + data_type = "landings", + uncertainty = get_data(data) |> dplyr::filter( .data$fleet == .env$fleets[i] & .data$type == "landings" ) |> - dplyr::select(-value) |> - dplyr::mutate( - estimation_type = "constant", - value = .data$uncertainty - ), - data_type = "landings" + dplyr::pull(dplyr::all_of("uncertainty")) ) } if ("age_comp" %in% fleet_types) { fleet_agecomp_distribution[[i]] <- initialize_data_distribution( module = fleet[[i]], - # TODO: need to update family and match options from the distribution - # column from the parameters tibble - family = multinomial(link = "logit"), - data_type = "agecomp" + data_type = "age_comp", + uncertainty = get_data(data) |> + dplyr::filter( + .data$fleet == .env$fleets[i] & + .data$type == "age_comp" + ) |> + dplyr::pull(dplyr::all_of("uncertainty")) ) } if ("length_comp" %in% fleet_types) { fleet_lengthcomp_distribution[[i]] <- initialize_data_distribution( module = fleet[[i]], - # TODO: need to update family and match options from the distribution - # column from the parameters tibble - family = multinomial(link = "logit"), - data_type = "lengthcomp" + data_type = "length_comp", + uncertainty = get_data(data) |> + dplyr::filter( + .data$fleet == .env$fleets[i] & + .data$type == "length_comp" + ) |> + dplyr::pull(dplyr::all_of("uncertainty")) ) } } diff --git a/R/parse_data_distribution.R b/R/parse_data_distribution.R new file mode 100644 index 000000000..638365a70 --- /dev/null +++ b/R/parse_data_distribution.R @@ -0,0 +1,168 @@ +#' Parse a distribution formula into structured components +#' +#' Parse a one-sided or two-sided formula-like distribution specification and +#' return the distribution family and parameter expressions as a list. This is +#' useful for translating user-written distribution declarations into a form +#' that can be validated or forwarded to downstream model-building utilities. +#' +#' @param expr A distribution specification provided as one of the following: +#' a character string (for example, `"~ gaussian(mean = x, sd = 0.1)"`), +#' an expression of length one (for example, +#' `expression(~ gaussian(mean = x, sd = 0.1))`), a language call produced +#' by `quote()`, or an R formula object. +#' +#' @details +#' The parser accepts inputs that represent a formula containing a distribution +#' call on the right-hand side, such as `~ gaussian(mean = x, sd = 0.1)`. +#' +#' Input handling and validation proceed in stages: +#' +#' * Character input is converted with [base::str2lang()]. +#' * `expression()` input must contain exactly one element, which is extracted +#' and parsed as the language object. +#' * Language calls and formula objects are used directly. +#' * Any other input type throws an error. +#' +#' After conversion, the object must be a call whose head is `~`. +#' Both one-sided formulas (`~ dist(...)`) and two-sided formulas +#' (`response ~ dist(...)`) are supported; in both cases, only the +#' right-hand side distribution call is parsed. +#' +#' The distribution call is decomposed into: +#' +#' * `family`: the function name used for the distribution, as a character +#' scalar (for example, `"gaussian"`, `"multinomial"`). +#' * `parameters`: a list of unevaluated parameter expressions from the +#' distribution call. Named arguments preserve their names, while positional +#' arguments remain unnamed. +#' +#' Parameter values are returned as language objects rather than evaluated +#' numeric vectors. This preserves references to symbols in the calling +#' environment (for example, `landings_expected`) so evaluation can occur later +#' in a context where those objects exist. +#' +#' @return A list with two elements: +#' +#' * `family`: character scalar containing the distribution family name. +#' * `parameters`: list of unevaluated arguments supplied to the distribution +#' call. +#' +#' @examples +#' FIMS:::parse_data_distribution( +#' "~ dmultinom(prob = age_comp_expected, size = 300)" +#' ) +#' FIMS:::parse_data_distribution( +#' quote(~ dnorm(mean = 0, sd = 1)) +#' ) +#' FIMS:::parse_data_distribution( +#' expression(~ dlnorm(meanlog = mu, sdlog = 0.1)) +#' ) +parse_data_distribution <- function(expr) { + # Preserve row alignment for upstream callers that keep NA placeholders. + if (is.atomic(expr) && length(expr) == 1L && is.na(expr)) { + all_parameter_names <- unique(unlist( + FIMS::fims_distributions$other_parameters + )) + all_parameters <- rep(1, length(all_parameter_names)) + names(all_parameters) <- all_parameter_names + return(c( + family = NA_character_, + as.list(all_parameters) + )) + } + + # Convert string to a language object if it isn't already + if (is.character(expr)) { + lang <- str2lang(expr) + } else if (is.expression(expr) && length(expr) == 1) { + lang <- expr[[1]] + } else if (is.call(expr) || inherits(expr, "formula")) { + lang <- expr + } else { + cli::cli_abort(".arg{expr} must be a string, expression, call, or formula.") + } + if (!is.call(lang)) { + cli::cli_abort("You must pass a string call.") + } + if (lang[[1]] != quote(`~`)) { + cli::cli_abort(c( + "Expression must be a distribution call", + "i" = "(e.g., ~dnorm(mean = landings_expected, sd = 1))" + )) + } + if (length(lang) == 3) { + cli::cli_abort(c( + "Two-sided formulas are not allowed in the data specifications", + "i" = "You passed {.code {lang}}" + )) + } + right_hand <- lang[[2]] + + # Extract the distribution family name + dist_name <- as.character(right_hand[[1]]) + dist_formals <- names(formals(dist_name)) + args_list <- as.list(right_hand[-1]) + link_name <- dplyr::filter( + FIMS::fims_distributions, + family == dist_name + ) |> + dplyr::pull(.data$central_parameter) + + # Checks + if (!dist_name %in% FIMS::fims_distributions[["family"]]) { + cli::cli_abort(c( + x = "The distribution you specified, {.code {dist_name}} is not available + in FIMS", + i = "Available distributions include + {.code {FIMS::fims_distributions$family}}." + )) + } + if (length(args_list[[link_name]]) != 1) { + cli::cli_abort(c( + "The central parameter does not exist for {dist_name} in your formula", + "You passed {expr}", + "The {dist_name} expects {.code {link_name}} as a parameter" + )) + } + if (length(setdiff(names(args_list), dist_formals)) > 0) { + cli::cli_abort(c( + x = "The arguments you passed in your distribution, {expr}, are not + formal arguments of the {.function {dist_name} function}", + i = "Argument options are {dist_formals}", + i = "Run {.code ?{dist_name}} for information on what each argument means" + )) + } + if (dist_name == "dnorm" && eval(args_list[["sd"]]) <= 0) { + cli::cli_abort(c( + x = "{.var sd} must be positive", + i = "In your distribution, {expr}, you passed {eval(args_list$sd)}" + )) + } + if (dist_name == "dlnorm" && eval(args_list[["sdlog"]]) <= 0) { + cli::cli_abort(c( + x = "{.var sdlog} must be positive", + i = "In your distribution, {expr}, you passed {eval(args_list$sdlog)}" + )) + } + + + return(c( + family = dist_name, + list(link = args_list[[link_name]]), + args_list + )) +} + +transpose_data_distribution <- function(x_list) { + data <- purrr::transpose(x_list) + unique_family <- unique(stats::na.omit(unlist(data[["family"]]))) + if (length(unique_family) > 1) { + cli::cli_abort(c( + "!" = "Cannot accommodate more than one family per data type right now.", + x = "You passed the following types: {unique_family}" + )) + } + + # TODO: Check that the mean/prob args have a single value + data +} diff --git a/R/setup_default_parameters.R b/R/setup_default_parameters.R index 8ac32eb50..332a48939 100644 --- a/R/setup_default_parameters.R +++ b/R/setup_default_parameters.R @@ -9,6 +9,12 @@ #' helper functions used within this function. The resulting tibble will have #' all of the necessary parameters, specific to your data, to run a FIMS #' model. +#' Initializing a FIMS model, i.e., `initialize_fims()` requires a `FIMSFrame` +#' object and a tibble of parameters. The parameter tibble can be automatically +#' generated using this function or by building up your own tibble from +#' helper functions used within this function. The resulting tibble will have +#' all of the necessary parameters, specific to your data, to run a FIMS +#' model. #' #' @details #' The function builds module-specific defaults by calling helper @@ -39,7 +45,7 @@ #' \item{\code{label}:}{The parameter name (e.g., "inflection_point").} #' \item{\code{age}:}{The age the parameter applies to.} #' \item{\code{length}:}{The length bin the parameter applies to.} -#' \item{\code{time}:}{The time step (year) the parameter applies to.} +#' \item{\code{timing}:}{The timing step (year) the parameter applies to.} #' \item{\code{value}:}{The initial value of the parameter.} #' \item{\code{estimation_type}:}{The estimation type (e.g., "constant", #' "fixed_effects", "random_effects").} @@ -243,7 +249,7 @@ setup_default_parameters_template <- function(n_parameters = 1) { label = NA_character_, age = NA_real_, length = NA_real_, - time = NA_integer_, + timing = NA_integer_, value = NA_real_, estimation_type = NA_character_, distribution_type = NA_character_, @@ -377,7 +383,7 @@ setup_default_Population <- function( label = "log_M", value = .env$log_M, age = rep(ages, n_years), - time = rep(years, each = n_ages), + timing = rep(years, each = n_ages), estimation_type = "constant" ) |> dplyr::add_row( @@ -568,7 +574,7 @@ setup_default_Fleet <- function( module_name = "Fleet", fleet = .env$fleet, label = "log_Fmort", - time = get_start_year(data):get_end_year(data), + timing = get_start_year(data):get_end_year(data), value = if (has_index) -200 else -3, estimation_type = if (has_index) "constant" else "fixed_effects" ) @@ -679,7 +685,7 @@ setup_default_BevertonHoltRecruitment <- function( # TODO: should this be log_recruit_dev to match output? label = "log_devs", value = 0.0, - time = (get_start_year(data) + 1):get_end_year(data), + timing = (get_start_year(data) + 1):get_end_year(data), estimation_type = "random_effects", distribution_type = "process", distribution = .env$distribution @@ -689,11 +695,11 @@ setup_default_BevertonHoltRecruitment <- function( dplyr::rows_update( tibble::tibble( label = "log_devs", - time = (get_start_year(data) + 1):get_end_year(data), + timing = (get_start_year(data) + 1):get_end_year(data), estimation_type = "constant", distribution_type = NA_character_ ), - by = c("label", "time") + by = c("label", "timing") ) } diff --git a/R/tidy-glance-fims.R b/R/tidy-glance-fims.R index c7a8abe10..11391ddb0 100644 --- a/R/tidy-glance-fims.R +++ b/R/tidy-glance-fims.R @@ -195,7 +195,7 @@ tidy.FIMSFit <- function( #' extracted from the TMB report. Returns a list-column when multiple #' populations are present.} #' \item{`fims_version`}{The version of FIMS used to fit the model.} -#' \item{`runtime_secs`}{Total wall-clock time of the fit in seconds.} +#' \item{`run_time`}{Total wall-clock time of the fit in seconds.} #' } #' #' @examples @@ -272,8 +272,8 @@ glance.FIMSFit <- function(x, ...) { } # run time - timing <- get_timing(x) - runtime_secs <- as.numeric(timing[["time_total"]], units = "secs") + run_time <- get_run_time(x) + run_time <- as.numeric(run_time[["time_total"]], units = "secs") tibble::tibble( logLik = log_lik, @@ -288,6 +288,6 @@ glance.FIMSFit <- function(x, ...) { converged = converged, terminal_ssb = terminal_ssb, fims_version = as.character(get_version(x)), - runtime_secs = runtime_secs + run_time = run_time ) } diff --git a/data-raw/data_big.R b/data-raw/data_big.R index 20f7ae332..a1311307f 100644 --- a/data-raw/data_big.R +++ b/data-raw/data_big.R @@ -261,7 +261,11 @@ landings_data <- data.frame( timing = returned_om[["om_input"]][["year"]], value = returned_om[["em_input"]][["L.obs"]][[1]], unit = "mt", # metric tons - uncertainty = cv_2_sd(returned_om[["em_input"]][["cv.L"]][[1]]) + uncertainty = paste( + "~ dlnorm(meanlog = log_landings_expected, sdlog = ", + cv_2_sd(returned_om[["em_input"]][["cv.L"]][[1]]), + ")" + ) ) ############################################################################### @@ -274,7 +278,11 @@ index_data <- data.frame( timing = returned_om[["om_input"]][["year"]], value = returned_om[["em_input"]][["surveyB.obs"]][[1]], unit = "mt", - uncertainty = cv_2_sd(returned_om[["em_input"]][["cv.survey"]][[1]]) + uncertainty = paste( + "~ dlnorm(meanlog = log_index_expected, sdlog = ", + cv_2_sd(returned_om[["em_input"]][["cv.survey"]][[1]]), + ")" + ) ) ############################################################################### @@ -285,14 +293,22 @@ age_data <- rbind( fleet = names(returned_om[["em_input"]][["n.L"]]), returned_om[["em_input"]][["L.age.obs"]][["fleet1"]], unit = "proportion", - uncertainty = returned_om[["em_input"]][["n.L"]][["fleet1"]], + uncertainty = paste( + "~ dmultinom(prob = agecomp_proportion, size = ", + returned_om[["em_input"]][["n.L"]][["fleet1"]], + ")" + ), timing = returned_om[["om_input"]][["year"]] ), data.frame( fleet = names(returned_om[["om_output"]][["survey_age_comp"]])[1], returned_om[["em_input"]][["survey.age.obs"]][[1]], unit = "proportion", - uncertainty = returned_om[["om_input"]][["n.survey"]][["survey1"]], + uncertainty = paste( + "~ dmultinom(prob = agecomp_proportion, size = ", + returned_om[["om_input"]][["n.survey"]][["survey1"]], + ")" + ), timing = returned_om[["om_input"]][["year"]] ) ) |> @@ -358,10 +374,7 @@ length_age_data <- data.frame( timing = NA_integer_, value = c(t(returned_om[["em_input"]][["age_to_length_conversion"]])), unit = "proportion", - uncertainty = c( - em_input[["n.L.lengthcomp"]][["fleet1"]], - em_input[["n.survey.lengthcomp"]][["survey1"]] - ) + uncertainty = NA_character_ ) # Create a length-composition data frame that will be filled by transforming @@ -377,12 +390,16 @@ length_comp_data <- data.frame( c(t(returned_om[["em_input"]][["survey.length.obs"]][["survey1"]])) ), unit = "proportion", - uncertainty = rep( - c( - em_input[["n.L.lengthcomp"]][["fleet1"]], - em_input[["n.survey.lengthcomp"]][["survey1"]] + uncertainty = paste( + "~ dmultinom(prob = lengthcomp_proportion, size = ", + rep( + c( + em_input[["n.L.lengthcomp"]][["fleet1"]], + em_input[["n.survey.lengthcomp"]][["survey1"]] + ), + length(len_bins) * length(timing_fishery[["timing"]]) ), - length(len_bins) * length(timing_fishery[["timing"]]) + ")" ) ) diff --git a/data-raw/fims_distributions.R b/data-raw/fims_distributions.R new file mode 100644 index 000000000..393fe9bb5 --- /dev/null +++ b/data-raw/fims_distributions.R @@ -0,0 +1,9 @@ +fims_distributions <- tibble::tribble( + ~family, ~sample_function, ~central_parameter, ~other_parameters, + "dnorm", "rnorm", "mean", list("sd"), + "dlnorm", "rlnorm", "meanlog", list("sdlog"), + "dmultinom", "rmultinom", "prob", list("size") +) + +usethis::use_data(fims_distributions, overwrite = TRUE) +on.exit(rm(fims_distributions), add = TRUE) diff --git a/data-raw/fims_input_types.R b/data-raw/fims_input_types.R index 68a544d54..ec223c8f9 100644 --- a/data-raw/fims_input_types.R +++ b/data-raw/fims_input_types.R @@ -1,5 +1,5 @@ fims_input_types <- c( - "age_comp", + "AgeComp", "age_to_length_conversion", "index", # TODO: change to singular landing diff --git a/data/data_big.rda b/data/data_big.rda index 6c58492a9..92641f5d1 100644 Binary files a/data/data_big.rda and b/data/data_big.rda differ diff --git a/data/fims_distributions.rda b/data/fims_distributions.rda new file mode 100644 index 000000000..897958586 Binary files /dev/null and b/data/fims_distributions.rda differ diff --git a/inst/include/README.md b/inst/include/README.md index 75bca885f..beafc9a81 100644 --- a/inst/include/README.md +++ b/inst/include/README.md @@ -42,7 +42,7 @@ The `TMB.h` file is the main wrapper around TMB's core headers. It supports both The `config.h` file provides the TMB/RTMB configuration used by FIMS. It enables Rcpp-compatible exception handling for TMB failures, turns on `TMB_SAFEBOUNDS`, selects the `TMBAD_FRAMEWORK`, uses a 64-bit index type (`uint64_t`), and configures thread-safe R support when OpenMP is available. -The `init_tmb.hpp` file handles runtime initialization and registration for the compiled shared object. It provides the `R_init_FIMS` registration callback, supports optional initialization through `FIMS_ONLOAD_INIT_TMB`, and registers TMB C-callables when `TMB_CCALLABLES` is available. +The `init_tmb.hpp` file handles run-time initialization and registration for the compiled shared object. It provides the `R_init_FIMS` registration callback, supports optional initialization through `FIMS_ONLOAD_INIT_TMB`, and registers TMB C-callables when `TMB_CCALLABLES` is available. ### inst/include/models/ diff --git a/man/FIMSFit.Rd b/man/FIMSFit.Rd index 3f5fdcdf6..c51855df2 100644 --- a/man/FIMSFit.Rd +++ b/man/FIMSFit.Rd @@ -9,7 +9,7 @@ FIMSFit( obj, opt = list(), sdreport = list(), - timing = c(time_total = as.difftime(0, units = "secs")), + run_time = c(time_total = as.difftime(0, units = "secs")), version = utils::packageVersion("FIMS") ) } @@ -24,7 +24,7 @@ FIMSFit( \item{sdreport}{An object of the \code{sdreport} class as returned from \code{\link[TMB:sdreport]{TMB::sdreport()}}.} -\item{timing}{A vector of at least length one, where all entries are of the +\item{run_time}{A vector of at least length one, where all entries are of the \code{timediff} class and at least one is named "time_total". This information is available in \code{\link[=fit_fims]{fit_fims()}} and added to this argument internally but if you are a power user you can calculate the time it took to run your model @@ -58,7 +58,7 @@ A list containing the model report from \code{obj[["report"]]()}. An object with the \code{sdreport} class containing the output from \code{TMB::sdreport(obj)}. } -\item{\code{timing}:}{ +\item{\code{run_time}:}{ The length of time it took to run the model if it was optimized. } \item{\code{version}:}{ diff --git a/man/FIMSFrame.Rd b/man/FIMSFrame.Rd index 6fbff87aa..2e7929abc 100644 --- a/man/FIMSFrame.Rd +++ b/man/FIMSFrame.Rd @@ -28,7 +28,8 @@ needed for different types of models. \subsection{data}{ The input data are both sorted (see the section below on sorting) and -expanded before returning them in the data slot. +expanded to include -999 values for all missing rows before returning them +in the data slot. \subsection{Ages}{ Currently, ages must be integers, i.e., FIMS cannot accommodate numeric ages @@ -43,6 +44,31 @@ multiple years, you must have age-2, age-3, and age-4 not just age-2 and age-4 fish in your composition data. } +\subsection{Uncertainty}{ + +Uncertainty information for your data contains information for fitting the +model and for creating bootstrapped data sets. Right-handed formulas are used +to specify the distribution your data is assumed to follow, the derived +quantity in the model that the data is assumed to represent, and parameters +associated with the assumed distribution. For example, landings data can be +distributed using a normal or lognormal distribution and if you use a +normal distribution the value is assumed to represent \code{landings_expected} but +if you use a lognormal distribution the value is assumed to represent +\code{log_landings_expected}. Therefore, for the former you would specify +\code{"~dnorm(mean = landings_expected, sd = 0.01)"} and for the later you would +specify \code{"~dlnorm(meanlog = log_landings_expected, sdlog = 0.01)"}. Where, +the input value for sd or sdlog is the standard deviation of the mean value +or standard deviation of the log of the mean value, respectively. The +distribution names and argument names match the names used in R. Use +\code{\link[=args]{args()}}, e.g., \code{args(dnorm)}, to get information on the argument names that +are required for your distribution. And, use \code{?distributions} to get a list +of distributions that have density functions in base R. To get information on +available derived quantities to set the central tendency of the distribution +to, you can create a fleet object and inspect the available names, e.g., +\code{names(methods::new(Fleet))}, or you can navigate to the Public Attributes +section of the doxygen documentation for Fleet by running \code{?Fleet}. +} + \subsection{Sorting}{ It is important that the order of the rows in the data are correct but it is diff --git a/man/data_big.Rd b/man/data_big.Rd index 2702a9621..9fdf24a9c 100644 --- a/man/data_big.Rd +++ b/man/data_big.Rd @@ -27,11 +27,9 @@ units for each data type are as follows. \code{mt} is used for \code{index}, \code{landings}, and \code{weight_at_age} data. \code{number} or \code{proportion} are each viable units for the composition data, where the former is the preferred unit of measurement.} -\item{uncertainty}{A real value providing a measurement of uncertainty -for value. For landings and indices of abundance this should be -the standard deviation of the logged observations if you are using the -lognormal distribution to fit your data. For composition data it will -be your input sample size. +\item{uncertainty}{A right-handed formula specifying the distributional +assumptions for the entry in \code{value}. See \code{\link[=FIMSFrame]{FIMSFrame()}} for more +information about this column. } } } @@ -49,4 +47,9 @@ The source code to make this dataset is in \code{data-raw/data_big.R} which is o GitHub but will not be in your local installation because it is in the .Rbuildignore file. } +\seealso{ +\itemize{ +\item \code{\link[=FIMSFrame]{FIMSFrame()}} +} +} \keyword{datasets} diff --git a/man/fims_distributions.Rd b/man/fims_distributions.Rd new file mode 100644 index 000000000..cae9f8aaa --- /dev/null +++ b/man/fims_distributions.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/fims_distributions.R +\docType{data} +\name{fims_distributions} +\alias{fims_distributions} +\title{Distribution parameter lookup table} +\format{ +A data frame with columns \code{family}, \code{sample_function}, \code{central_parameter}, and +\code{other_parameters}. \code{sample_function} is a character column with function +names. +} +\usage{ +fims_distributions +} +\description{ +A lookup table that maps each supported distribution family to the field +name used to link observed data, the random-value generator, and the +remaining parameter names for that family. +} +\keyword{datasets} diff --git a/man/get_FIMSFit.Rd b/man/get_FIMSFit.Rd index 439aa4405..54eb2225b 100644 --- a/man/get_FIMSFit.Rd +++ b/man/get_FIMSFit.Rd @@ -20,8 +20,8 @@ \alias{get_estimates,FIMSFit-method} \alias{get_number_of_parameters} \alias{get_number_of_parameters,FIMSFit-method} -\alias{get_timing} -\alias{get_timing,FIMSFit-method} +\alias{get_run_time} +\alias{get_run_time,FIMSFit-method} \alias{get_version} \alias{get_version,FIMSFit-method} \alias{get_model_output} @@ -64,9 +64,9 @@ get_number_of_parameters(x) \S4method{get_number_of_parameters}{FIMSFit}(x) -get_timing(x) +get_run_time(x) -\S4method{get_timing}{FIMSFit}(x) +\S4method{get_run_time}{FIMSFit}(x) get_version(x) @@ -108,7 +108,7 @@ uncertainties from a fitted model. number of fixed-effect parameters and the number of random-effect parameters in the model. -\code{\link[=get_timing]{get_timing()}} returns the amount of time it took to run the model in +\code{\link[=get_run_time]{get_run_time()}} returns the amount of time it took to run the model in seconds as a \code{difftime} object. \code{\link[=get_version]{get_version()}} returns the \code{package_version} of FIMS that was used to fit diff --git a/man/glance.FIMSFit.Rd b/man/glance.FIMSFit.Rd index 2ce8235e4..b249a3cfc 100644 --- a/man/glance.FIMSFit.Rd +++ b/man/glance.FIMSFit.Rd @@ -34,7 +34,7 @@ below 0.001 are generally considered well-converged.} extracted from the TMB report. Returns a list-column when multiple populations are present.} \item{\code{fims_version}}{The version of FIMS used to fit the model.} -\item{\code{runtime_secs}}{Total wall-clock time of the fit in seconds.} +\item{\code{run_time}}{Total wall-clock time of the fit in seconds.} } } \description{ diff --git a/man/initialize_data_distribution.Rd b/man/initialize_data_distribution.Rd index 2e1218199..b53829f89 100644 --- a/man/initialize_data_distribution.Rd +++ b/man/initialize_data_distribution.Rd @@ -8,9 +8,8 @@ \usage{ initialize_data_distribution( module, - family = NULL, - sd = tibble::tibble(value = 1, estimation_type = "constant"), - data_type = c("landings", "index", "agecomp", "lengthcomp") + data_type = c("landings", "index", "age_comp", "length_comp"), + uncertainty ) initialize_process_distribution( @@ -26,29 +25,15 @@ initialize_process_structure(module, par) \item{module}{An identifier to a C++ fleet module that is linked to the data of interest.} -\item{family}{A description of the error distribution and link function to -be used in the model. The argument takes a family class, e.g., -\code{stats::gaussian(link = "identity")}.} - -\item{sd}{A list of length two. The first entry is named \code{"value"} and it -stores the initial values (scalar or vector) for the relevant standard -deviations. The default is \code{value = 1}. The second entry is named -\code{"estimation_type"} and it stores a vector of booleans (default = -"constant") is a string indicating whether or not standard deviation is -estimated as a fixed effect or held constant. If \code{"value"} is a vector and -\code{"estimation_type"} is a scalar, the single value specified -\code{"estimation_type"} value will be repeated to match the length of \code{value}. -Otherwise, the dimensions of the two must match.} - \item{data_type}{A string specifying the type of data that the distribution will be fit to. Allowable types include -"landings", "index", "agecomp", and "lengthcomp" +"landings", "index", "age_comp", and "length_comp" and the default is landings.} -\item{par}{A string specifying the parameter name the distribution applies -to. Parameters must be members of the specified module. Use -\code{methods::show(module)} to obtain names of parameters within the module.} +\item{uncertainty}{A vector of strings specifying formulas for each data +point. See \code{\link[=FIMSFrame]{FIMSFrame()}} for more information on what the formula should +look like.} } \value{ A reference class. is returned. Use \code{\link[methods:show]{methods::show()}} to view the various @@ -70,12 +55,11 @@ n_years <- 30 # Create a new fleet module fleet <- methods::new(Fleet) # Create a distribution for the fleet module +sd_log <- rep(sqrt(log(0.01^2 + 1)), n_years) fleet_distribution <- initialize_data_distribution( module = fishing_fleet, - family = lognormal(link = "log"), - sd = list( - value = rep(sqrt(log(0.01^2 + 1)), n_years), - estimation_type = rep("constant", n_years) # Can be a single "constant" + uncertainty = glue::glue( + "~dlnorm(meanlog = log_index_expected, sdlog = {sd_log})" ), data_type = "index" ) diff --git a/man/parse_data_distribution.Rd b/man/parse_data_distribution.Rd new file mode 100644 index 000000000..460d4ba85 --- /dev/null +++ b/man/parse_data_distribution.Rd @@ -0,0 +1,72 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parse_data_distribution.R +\name{parse_data_distribution} +\alias{parse_data_distribution} +\title{Parse a distribution formula into structured components} +\usage{ +parse_data_distribution(expr) +} +\arguments{ +\item{expr}{A distribution specification provided as one of the following: +a character string (for example, \code{"~ gaussian(mean = x, sd = 0.1)"}), +an expression of length one (for example, +\code{expression(~ gaussian(mean = x, sd = 0.1))}), a language call produced +by \code{quote()}, or an R formula object.} +} +\value{ +A list with two elements: +\itemize{ +\item \code{family}: character scalar containing the distribution family name. +\item \code{parameters}: list of unevaluated arguments supplied to the distribution +call. +} +} +\description{ +Parse a one-sided or two-sided formula-like distribution specification and +return the distribution family and parameter expressions as a list. This is +useful for translating user-written distribution declarations into a form +that can be validated or forwarded to downstream model-building utilities. +} +\details{ +The parser accepts inputs that represent a formula containing a distribution +call on the right-hand side, such as \code{~ gaussian(mean = x, sd = 0.1)}. + +Input handling and validation proceed in stages: +\itemize{ +\item Character input is converted with \code{\link[base:str2lang]{base::str2lang()}}. +\item \code{expression()} input must contain exactly one element, which is extracted +and parsed as the language object. +\item Language calls and formula objects are used directly. +\item Any other input type throws an error. +} + +After conversion, the object must be a call whose head is \code{~}. +Both one-sided formulas (\code{~ dist(...)}) and two-sided formulas +(\code{response ~ dist(...)}) are supported; in both cases, only the +right-hand side distribution call is parsed. + +The distribution call is decomposed into: +\itemize{ +\item \code{family}: the function name used for the distribution, as a character +scalar (for example, \code{"gaussian"}, \code{"multinomial"}). +\item \code{parameters}: a list of unevaluated parameter expressions from the +distribution call. Named arguments preserve their names, while positional +arguments remain unnamed. +} + +Parameter values are returned as language objects rather than evaluated +numeric vectors. This preserves references to symbols in the calling +environment (for example, \code{landings_expected}) so evaluation can occur later +in a context where those objects exist. +} +\examples{ +FIMS:::parse_data_distribution( + "~ dmultinom(prob = age_comp_expected, size = 300)" +) +FIMS:::parse_data_distribution( + quote(~ dnorm(mean = 0, sd = 1)) +) +FIMS:::parse_data_distribution( + expression(~ dlnorm(meanlog = mu, sdlog = 0.1)) +) +} diff --git a/man/setup_default_parameters.Rd b/man/setup_default_parameters.Rd index 4827fd983..28fcabcc4 100644 --- a/man/setup_default_parameters.Rd +++ b/man/setup_default_parameters.Rd @@ -26,7 +26,7 @@ in PascalCase to match the names used in the C++ code.} \item{\code{label}:}{The parameter name (e.g., "inflection_point").} \item{\code{age}:}{The age the parameter applies to.} \item{\code{length}:}{The length bin the parameter applies to.} -\item{\code{time}:}{The time step (year) the parameter applies to.} +\item{\code{timing}:}{The timing step (year) the parameter applies to.} \item{\code{value}:}{The initial value of the parameter.} \item{\code{estimation_type}:}{The estimation type (e.g., "constant", "fixed_effects", "random_effects").} @@ -45,6 +45,12 @@ generated using this function or by building up your own tibble from helper functions used within this function. The resulting tibble will have all of the necessary parameters, specific to your data, to run a FIMS model. +Initializing a FIMS model, i.e., \code{initialize_fims()} requires a \code{FIMSFrame} +object and a tibble of parameters. The parameter tibble can be automatically +generated using this function or by building up your own tibble from +helper functions used within this function. The resulting tibble will have +all of the necessary parameters, specific to your data, to run a FIMS +model. } \details{ The function builds module-specific defaults by calling helper diff --git a/pkgdown/assets/fims-cheatsheet-thumb.png b/pkgdown/assets/fims-cheatsheet-thumb.png index 6ccd859c4..306f36b41 100644 Binary files a/pkgdown/assets/fims-cheatsheet-thumb.png and b/pkgdown/assets/fims-cheatsheet-thumb.png differ diff --git a/pkgdown/assets/fims-cheatsheet.pdf b/pkgdown/assets/fims-cheatsheet.pdf index 5435589b5..21af8be45 100644 Binary files a/pkgdown/assets/fims-cheatsheet.pdf and b/pkgdown/assets/fims-cheatsheet.pdf differ diff --git a/tests/testthat/_snaps/setup_default_parameters/default_parameters.csv b/tests/testthat/_snaps/setup_default_parameters/default_parameters.csv index 27f5ffee8..c544016dc 100644 --- a/tests/testthat/_snaps/setup_default_parameters/default_parameters.csv +++ b/tests/testthat/_snaps/setup_default_parameters/default_parameters.csv @@ -1,4 +1,4 @@ -"","module_name","fleet","module_type","label","age","length","time","value","estimation_type","distribution_type","distribution" +"","module_name","fleet","module_type","label","age","length","timing","value","estimation_type","distribution_type","distribution" "1","Selectivity","fleet1","Logistic","inflection_point",NA,NA,NA,2,"fixed_effects",NA,NA "2","Selectivity","fleet1","Logistic","slope",NA,NA,NA,1,"fixed_effects",NA,NA "3","Selectivity","survey1","Logistic","inflection_point",NA,NA,NA,2,"fixed_effects",NA,NA diff --git a/tests/testthat/fixtures/integration_test_data.RData b/tests/testthat/fixtures/integration_test_data.RData index 974e99181..ae6799bba 100644 Binary files a/tests/testthat/fixtures/integration_test_data.RData and b/tests/testthat/fixtures/integration_test_data.RData differ diff --git a/tests/testthat/fixtures/integration_test_data_components.RData b/tests/testthat/fixtures/integration_test_data_components.RData index ff6df14d6..f2efa2085 100644 Binary files a/tests/testthat/fixtures/integration_test_data_components.RData and b/tests/testthat/fixtures/integration_test_data_components.RData differ diff --git a/tests/testthat/helper-integration-tests-setup-run.R b/tests/testthat/helper-integration-tests-setup-run.R index 58a4fff46..5d9cc9f20 100644 --- a/tests/testthat/helper-integration-tests-setup-run.R +++ b/tests/testthat/helper-integration-tests-setup-run.R @@ -133,10 +133,10 @@ prepare_test_data <- function() { tibble::tibble( fleet = "fleet1", label = "log_Fmort", - time = 1:get_n_years(data_age_length_comp), + timing = 1:get_n_years(data_age_length_comp), value = log(om_output_list[[iter_id]][["f"]]), ), - by = c("fleet", "label", "time") + by = c("fleet", "label", "timing") ) |> # Update selectivity parameters and log_q for survey1 dplyr::rows_update( @@ -147,14 +147,14 @@ prepare_test_data <- function() { ), by = c("fleet", "label") ) |> - # Update log_devs in the Recruitment module (time steps 2-30) + # Update log_devs in the Recruitment module (timing steps 2-30) dplyr::rows_update( tibble::tibble( label = "log_devs", - time = 2:get_n_years(data_age_length_comp), + timing = 2:get_n_years(data_age_length_comp), value = om_input_list[[iter_id]][["logR.resid"]][-1] ), - by = c("label", "time") + by = c("label", "timing") ) |> # Update log_sd for log_devs in the Recruitment module # Note: logR_sd is the standard deviation on the natural scale of the diff --git a/tests/testthat/test-check_mle_convergence.R b/tests/testthat/test-check_mle_convergence.R index e74298e7e..2a0615d34 100644 --- a/tests/testthat/test-check_mle_convergence.R +++ b/tests/testthat/test-check_mle_convergence.R @@ -25,14 +25,14 @@ make_mock_mle_opt <- function(convergence = 0, message = NULL, par = numeric()) ) } -mock_fimsfit <- function(input, obj, opt, sdreport, timing) { +mock_fimsfit <- function(input, obj, opt, sdreport, run_time) { list( tag = "mock_fit", input = input, obj = obj, opt = opt, sdreport = sdreport, - timing = timing + run_time = run_time ) } diff --git a/tests/testthat/test-check_sdreport_convergence.R b/tests/testthat/test-check_sdreport_convergence.R index 9ca3155bc..dc14b83c9 100644 --- a/tests/testthat/test-check_sdreport_convergence.R +++ b/tests/testthat/test-check_sdreport_convergence.R @@ -122,14 +122,14 @@ test_that("check_sdreport_convergence() works with correct inputs", { out <- NULL testthat::with_mocked_bindings( - FIMSFit = function(input, obj, opt, sdreport, timing) { + FIMSFit = function(input, obj, opt, sdreport, run_time) { list( tag = "mock_fit", input = input, obj = obj, opt = opt, sdreport = sdreport, - timing = timing + run_time = run_time ) }, print = function(x, ...) invisible(x), diff --git a/tests/testthat/test-distribution-formulas.R b/tests/testthat/test-distribution-formulas.R index 39020819d..fe49c950c 100644 --- a/tests/testthat/test-distribution-formulas.R +++ b/tests/testthat/test-distribution-formulas.R @@ -73,14 +73,9 @@ fishing_fleet$SetObservedIndexDataID(fishing_fleet_index$get_id()) fleet_sd <- rep(sqrt(log(em_input$cv.L$fleet1^2 + 1)), om_input$nyr) fishing_fleet_index_distribution1 <- initialize_data_distribution( module = fishing_fleet, - family = lognormal(link = "log"), - sd = list(value = fleet_sd, estimation_type = "constant"), - data_type = "index" -) -fishing_fleet_index_distribution2 <- initialize_data_distribution( - module = fishing_fleet, - family = stats::gaussian(link = "log"), - sd = list(value = fleet_sd[1], estimation_type = "fixed_effects"), + uncertainty = glue::glue( + "~dlnorm(meanlog = log_index_expected, sdlog = {fleet_sd})" + ), data_type = "index" ) @@ -114,7 +109,11 @@ test_that("`initialize_data_distribution()` works with correct inputs", { #' @description Test that `initialize_data_distribution()` returns the correct log sd values when scalar. expect_equal( log(fleet_sd[1]), - fishing_fleet_index_distribution2$log_sd[1]$value + initialize_data_distribution( + fishing_fleet, + data_type = "index", + glue::glue("~dlnorm(meanlog = log_index_expected, sdlog = {fleet_sd[1]})") + )$log_sd[1]$value ) }) @@ -199,90 +198,23 @@ test_that("`initialize_process_distribution()` returns correct error messages", test_that("`initialize_data_distribution()` returns correct error messages", { - #' @description Test that error is thrown when `family` and `index` `data_type` don't match in `initialize_data_distribution()`. - expect_error( - initialize_data_distribution( - module = fishing_fleet, - family = multinomial(), - sd = list(value = fleet_sd, estimation_type = "constant"), - data_type = "index" - ), - "does not allow the family to be" - ) - - #' @description Test that error is thrown when `family` and `landings` `data_type` don't match in `initialize_data_distribution()`. - expect_error( - initialize_data_distribution( - module = fishing_fleet, - family = multinomial(), - sd = list(value = fleet_sd, estimation_type = "constant"), - data_type = "landings" - ), - "does not allow the family to be" - ) - - #' @description Test that error is thrown when `family` and `agecomp` `data_type` don't match in `initialize_data_distribution()`. - expect_error( - initialize_data_distribution( - module = fishing_fleet, - family = gaussian(), - sd = list(value = fleet_sd, estimation_type = "constant"), - data_type = "agecomp" - ), - "does not allow the family to be" - ) - - #' @description Test that error is thrown when `family` and `lengthcomp` `data_type` don't match in `initialize_data_distribution()`. - expect_error( - initialize_data_distribution( - module = fishing_fleet, - family = lognormal(), - sd = list(value = fleet_sd, estimation_type = "constant"), - data_type = "lengthcomp" - ), - "does not allow the family to be" - ) - #' @description Test that error is thrown when `data_type` is incorrect in `initialize_data_distribution()`. expect_error( initialize_data_distribution( module = fishing_fleet, - family = lognormal(), - sd = list(value = fleet_sd, estimation_type = "constant"), - data_type = "length_comp" + uncertainty = "~dmultinom(prob = agecomp_proportion, size = 10)", + data_type = "bad_type" ), "must be one of" ) - #' @description Test that error is thrown when `sd` value and `estimation_type` dimensions do not match. - expect_error( - initialize_data_distribution( - module = fishing_fleet, - family = multinomial(), - sd = list(value = fleet_sd, estimation_type = c("constant", "constant")), - data_type = "agecomp" - ), - "must match if more than one value" - ) - - #' @description Test that error is thrown when `sd` value is missing. - expect_error( - initialize_data_distribution( - module = fishing_fleet, - family = multinomial(), - sd = list(estimation_type = "constant"), - data_type = "agecomp" - ), - "need to be present" - ) - #' @description Test that error is thrown when `family` is missing. + #' @description Test that error is thrown when `uncertainty` is missing. expect_error( initialize_data_distribution( module = fishing_fleet, - sd = list(value = fleet_sd), - data_type = "agecomp" + data_type = "age_comp" ), - "is missing from" + "argument .*uncertainty.* is missing" ) clear() }) diff --git a/tests/testthat/test-fimsfit.R b/tests/testthat/test-fimsfit.R index effa49cc2..54e52891c 100644 --- a/tests/testthat/test-fimsfit.R +++ b/tests/testthat/test-fimsfit.R @@ -32,7 +32,7 @@ test_that("`is.FIMSFit()` works with correct inputs", { expected_names <- c( "input", "obj", "opt", "max_gradient", "gradient", "report", "sdreport", - "number_of_parameters", "timing", "version", "model_output" + "number_of_parameters", "run_time", "version", "model_output" ) #' @description Test a FIMSFit object has the correct slot names. expect_equal( @@ -46,19 +46,19 @@ test_that("`is.FIMSFit()` returns correct outputs for edge cases", { #' @description Test that `is.FIMSFit("not_a_FIMSFit")` returns FALSE. expect_false(is.FIMSFit("not_a_FIMSFit")) - # Modify the total time to be more than a day - fit_age_length_comp@timing[["time_total"]] <- 86401 # 60*60*24+1 - #' @description Test that `print(FIMSFit)` returns no error when the total time is more than a day. + # Modify the total run_time to be more than a day + fit_age_length_comp@run_time[["time_total"]] <- 86401 # 60*60*24+1 + #' @description Test that `print(FIMSFit)` returns no error when the total run_time is more than a day. expect_no_error(print(fit_age_length_comp)) - # Modify the total time to be more than a hour - fit_age_length_comp@timing[["time_total"]] <- 3601 # 60*60+1 - #' @description Test that `print(FIMSFit)` returns no error when the total time is more than an hour. + # Modify the total run_time to be more than a hour + fit_age_length_comp@run_time[["time_total"]] <- 3601 # 60*60+1 + #' @description Test that `print(FIMSFit)` returns no error when the total run_time is more than an hour. expect_no_error(print(fit_age_length_comp)) - # Modify the total time to be more than a minute - fit_age_length_comp@timing[["time_total"]] <- 61 # 60+1 - #' @description Test that `print(FIMSFit)` returns no error when the total time is more than a minute. + # Modify the total run_time to be more than a minute + fit_age_length_comp@run_time[["time_total"]] <- 61 # 60+1 + #' @description Test that `print(FIMSFit)` returns no error when the total run_time is more than a minute. expect_no_error(print(fit_age_length_comp)) }) @@ -131,10 +131,10 @@ test_that("fit_fims() errors when optimization fails to converge", { dplyr::rows_update( tibble::tibble( label = "log_devs", - time = 2:get_n_years(data_age_comp), + timing = 2:get_n_years(data_age_comp), estimation_type = "fixed_effects" ), - by = c("label", "time") + by = c("label", "timing") ) |> dplyr::rows_update( tibble::tibble( @@ -151,7 +151,7 @@ test_that("fit_fims() errors when optimization fails to converge", { tibble::tibble( fleet = "fleet1", type = "landings", - uncertainty = 10, + uncertainty = "~dlnorm(meanlog=log_landings_expected, sdlog = 10)", ), by = c("fleet", "type") ) |> @@ -170,10 +170,10 @@ test_that("fit_fims() errors when optimization fails to converge", { dplyr::rows_update( tibble::tibble( label = "log_devs", - time = 2:get_n_years(data_age_comp), + timing = 2:get_n_years(data_age_comp), estimation_type = "fixed_effects" ), - by = c("label", "time") + by = c("label", "timing") ) |> dplyr::rows_update( tibble::tibble( @@ -186,10 +186,10 @@ test_that("fit_fims() errors when optimization fails to converge", { dplyr::rows_update( tibble::tibble( label = "log_M", - time = 1, + timing = 1, estimation_type = "fixed_effects" ), - by = c("label", "time") + by = c("label", "timing") ) initialized_model <- parameters_4_model |> diff --git a/tests/testthat/test-get_timing.R b/tests/testthat/test-get_run_time.R similarity index 62% rename from tests/testthat/test-get_timing.R rename to tests/testthat/test-get_run_time.R index f7484cb0c..7f18bc886 100644 --- a/tests/testthat/test-get_timing.R +++ b/tests/testthat/test-get_run_time.R @@ -7,14 +7,14 @@ #' one lines, that will be used in the bookdown report of the results from #' {testthat}. This line can be more than 80 characters. -# get_timing ---- +# get_run_time ---- ## Setup ---- # Load or prepare any necessary data for testing if (!file.exists(testthat::test_path("fixtures", "fit_age_length_comp.RDS"))) { prepare_test_data() } ## IO correctness ---- -test_that("`get_timing()` works with correct inputs", { +test_that("`get_run_time()` works with correct inputs", { # Load the test data from an RDS file containing model fits. # List all RDS files in the fixtures directory that match the pattern "fit*_.RDS" fit_files <- list.files( @@ -28,32 +28,32 @@ test_that("`get_timing()` works with correct inputs", { ) # Function to read the RDS file and get input - check_timing <- function(fit_file) { + check_run_time <- function(fit_file) { fit_data <- readRDS(fit_file) - timing <- get_timing(fit_data) - #' @description Test that `get_timing()` returns correct output for the `timing` slot. + run_time <- get_run_time(fit_data) + #' @description Test that `get_run_time()` returns correct output for the `run_time` slot. expect_equal( - object = timing, - expected = fit_data@timing + object = run_time, + expected = fit_data@run_time ) - #' @description Test that `get_timing()` returns correct names for the `timing` slot. + #' @description Test that `get_run_time()` returns correct names for the `run_time` slot. expect_equal( - object = names(timing), + object = names(run_time), expected = expected_names ) - #' @description Test that `get_timing()` returns > 0 values for the `timing` slot. - expect_true(object = all(timing > 0)) + #' @description Test that `get_run_time()` returns > 0 values for the `run_time` slot. + expect_true(object = all(run_time > 0)) } # Use purrr::map to apply the function to each file - result <- purrr::map(fit_files, check_timing) + result <- purrr::map(fit_files, check_run_time) }) ## Edge handling ---- -test_that("`get_timing()` returns correct outputs for edge cases", { - #' @description Test that `get_timing()` returns an error when given invalid input. +test_that("`get_run_time()` returns correct outputs for edge cases", { + #' @description Test that `get_run_time()` returns an error when given invalid input. expect_error( - object = get_timing("invalid_input") + object = get_run_time("invalid_input") ) }) diff --git a/tests/testthat/test-initialize_modules.R b/tests/testthat/test-initialize_modules.R index e97779f9b..8872b85c8 100644 --- a/tests/testthat/test-initialize_modules.R +++ b/tests/testthat/test-initialize_modules.R @@ -41,17 +41,17 @@ test_that("`initialize_fims()` works with edge cases", { # to constant dplyr::mutate( estimation_type = dplyr::if_else( - time <= 11, + timing <= 11, "constant", estimation_type ), distribution_type = dplyr::if_else( - time <= 11, + timing <= 11, NA_character_, distribution_type ), distribution = dplyr::if_else( - time <= 11, + timing <= 11, NA_character_, distribution ) @@ -193,10 +193,10 @@ test_that("`initialize_fims()` returns correct error messages", { y = tibble::tibble( module_name = "Recruitment", label = "log_devs", - time = 2:get_n_years(data), + timing = 2:get_n_years(data), estimation_type = "random_effects" ), - by = c("module_name", "label", "time") + by = c("module_name", "label", "timing") ) #' @description Test that `initialize_fims()` returns correct error with distribution estimation type mismatch. @@ -213,10 +213,10 @@ test_that("`initialize_fims()` returns correct error messages", { y = tibble::tibble( module_name = "Recruitment", label = "log_devs", - time = 2:get_n_years(data), + timing = 2:get_n_years(data), estimation_type = "constant" ), - by = c("module_name", "label", "time") + by = c("module_name", "label", "timing") ) #' @description Test that `initialize_recruitment()` handles missing distribution for recruitment correctly. @@ -546,7 +546,13 @@ test_that("`initialize_comp()` works with correct inputs", { data |> get_data() |> dplyr::filter(type == "age_comp", fleet == "fleet1") |> - dplyr::mutate(out = value * uncertainty) |> + dplyr::mutate( + out = value * purrr::map_dbl( + purrr::map(.f = parse_data_distribution, .x = .data$uncertainty), + list("size"), + .default = NA_real_ + ) + ) |> dplyr::pull(out) ) clear() @@ -573,7 +579,11 @@ test_that("`initialize_comp()` works with correct inputs", { data |> get_data() |> dplyr::filter(type == "length_comp", fleet == "fleet1") |> - dplyr::mutate(out = value * uncertainty) |> + dplyr::mutate(out = value * purrr::map_dbl( + purrr::map(.f = parse_data_distribution, .x = .data$uncertainty), + list("size"), + .default = NA_real_ + )) |> dplyr::pull(out) ) clear() @@ -587,17 +597,17 @@ test_that("`initialize_fims()` works with edge cases", { # to constant dplyr::mutate( estimation_type = dplyr::if_else( - time <= 11, + timing <= 11, "constant", estimation_type ), distribution_type = dplyr::if_else( - time <= 11, + timing <= 11, NA_character_, distribution_type ), distribution = dplyr::if_else( - time <= 11, + timing <= 11, NA_character_, distribution ) diff --git a/tests/testthat/test-integration-caa-mle-with-wrappers.R b/tests/testthat/test-integration-caa-mle-with-wrappers.R index d964576a9..e2a556d32 100644 --- a/tests/testthat/test-integration-caa-mle-with-wrappers.R +++ b/tests/testthat/test-integration-caa-mle-with-wrappers.R @@ -220,7 +220,7 @@ test_that("catch-at-age model (estimation MLE with wrappers) works with mixed es fit_mixed_estimation_types <- modified_parameters |> dplyr::mutate( estimation_type = dplyr::if_else( - fleet == "fleet1" & label == "log_Fmort" & time %in% 1:10, + fleet == "fleet1" & label == "log_Fmort" & timing %in% 1:10, "constant", estimation_type ) @@ -267,11 +267,11 @@ test_that("catch-at-age model (estimation MLE with wrappers) returns an error wh # log_devs has a special error when set to constant y = tibble::tibble( label = "log_devs", - time = 2:get_n_years(data_age_length_comp), + timing = 2:get_n_years(data_age_length_comp), distribution_type = NA_character_, distribution = NA_character_, ), - by = c("label", "time") + by = c("label", "timing") ) |> dplyr::rows_update( # log_sd has a special error when there isn't a log_devs or log_r parameter set @@ -335,10 +335,10 @@ test_that("catch-at-age model (estimation MLE with wrappers) returns an error wh dplyr::rows_update( y = tibble::tibble( label = "log_devs", - time = 2:get_n_years(data_age_length_comp), + timing = 2:get_n_years(data_age_length_comp), estimation_type = "constant" ), - by = c("label", "time") + by = c("label", "timing") ) #' @description Test that FIMS returns an error when log_devs are constant but Recruitment expects a distribution process. diff --git a/tests/testthat/test-integration-fims-estimation-random-effects-with-wrappers.R b/tests/testthat/test-integration-fims-estimation-random-effects-with-wrappers.R index 4a562ac78..fcba78783 100644 --- a/tests/testthat/test-integration-fims-estimation-random-effects-with-wrappers.R +++ b/tests/testthat/test-integration-fims-estimation-random-effects-with-wrappers.R @@ -502,9 +502,9 @@ test_that("estimation test with recruitment re on logr", { #' @description Test that the `expected_recruitment` from both `fit_log_r` and `fit_log_devs` runs are approximately equal within a tolerance of 0.001. expect_equal(fit_log_r@report[["expected_recruitment"]], fit_log_devs@report[["expected_recruitment"]], tolerance = .001) #' @description Test that the `time_optimization` from `fit_log_r` is less than or equal to that from `fit_log_devs`. - expect_lte(fit_log_r@timing[["time_optimization"]], fit_log_devs@timing[["time_optimization"]]) + expect_lte(fit_log_r@run_time[["time_optimization"]], fit_log_devs@run_time[["time_optimization"]]) #' @description Test that the `time_sdreport` from `fit_log_r` is less than or equal to that from `fit_log_devs`. - expect_lte(fit_log_r@timing[["time_sdreport"]], fit_log_devs@timing[["time_sdreport"]]) + expect_lte(fit_log_r@run_time[["time_sdreport"]], fit_log_devs@run_time[["time_sdreport"]]) clear() }) diff --git a/tests/testthat/test-parse_data_distribution.R b/tests/testthat/test-parse_data_distribution.R new file mode 100644 index 000000000..f575e2179 --- /dev/null +++ b/tests/testthat/test-parse_data_distribution.R @@ -0,0 +1,106 @@ +# Instructions ---- +#' This file follows the format generated by FIMS:::use_testthat_template(). +#' Necessary tests include input and output (IO) correctness [IO +#' correctness], edge-case handling [Edge handling], and built-in errors and +#' warnings [Error handling]. See `?FIMS:::use_testthat_template` for more +#' information. Every test should have a @description tag, which can only span +#' one lines, that will be used in the bookdown report of the results from +#' {testthat}. This line can be more than 80 characters. + +# parse_data_distribution ---- +## Setup ---- +# No setup required. + +## IO correctness ---- +test_that("`parse_data_distribution()` works with correct inputs", { + #' @description Test that a lognormal distribution string is parsed into a flat list of its family, link, and parameters. + expect_equal( + parse_data_distribution( + "~ dlnorm(meanlog = log_landings_expected, sdlog = 0.00999975001354021)" + ), + list( + family = "dlnorm", + link = quote(log_landings_expected), + meanlog = quote(log_landings_expected), + sdlog = 0.00999975001354021 + ) + ) + + #' @description Test that whitespace is not required between the formula operator and a normal distribution call. + expect_equal( + parse_data_distribution("~dnorm(mean = landings_expected, sd = 0.1)"), + list( + family = "dnorm", + link = quote(landings_expected), + mean = quote(landings_expected), + sd = 0.1 + ) + ) + + #' @description Test that calls and single-element expressions are accepted without first converting them to strings. + expect_equal( + parse_data_distribution(quote(~dnorm(mean = index_expected, sd = 1))), + parse_data_distribution(expression(~dnorm(mean = index_expected, sd = 1))) + ) +}) + +## Edge handling ---- +test_that("`parse_data_distribution()` returns correct outputs for edge cases", { + #' @description Test that a missing distribution specification returns an `NA` family while preserving row alignment. + expect_equal( + parse_data_distribution(NA_character_)[1], + list(family = NA_character_) + ) +}) + +## Error handling ---- +test_that("`parse_data_distribution()` returns correct error messages", { + #' @description Test that unsupported input types return an informative error. + expect_error( + parse_data_distribution(list("~dnorm(mean = x, sd = 1)")), + "must be a string, expression, call, or formula", + fixed = TRUE + ) + + #' @description Test that input without the formula operator returns an informative error. + expect_error( + parse_data_distribution("dnorm(mean = x, sd = 1)"), + "Expression must be a distribution call", + fixed = TRUE + ) + + #' @description Test that two-sided formulas are rejected for data distribution specifications. + expect_error( + parse_data_distribution("landings ~ dnorm(mean = x, sd = 1)"), + "Two-sided formulas are not allowed", + fixed = TRUE + ) + + #' @description Test that the family-specific central parameter must be present. + expect_error( + parse_data_distribution("~dnorm(sd = 1)"), + "The central parameter does not exist for dnorm", + fixed = TRUE + ) + + #' @description Test that distribution arguments must be formal arguments of the selected family. + expect_error( + parse_data_distribution("~dnorm(mean = x, sd = 1, extra = 2)"), + "are not formal arguments of the dnorm function", + fixed = TRUE + ) + + #' @description Test that a normal distribution standard deviation must be positive. + expect_error( + parse_data_distribution("~dnorm(mean = x, sd = -0.1)"), + "`sd` must be positive", + fixed = TRUE + ) + + #' @description Test that a lognormal distribution log-scale standard deviation must be positive. + expect_error( + parse_data_distribution("~dlnorm(meanlog = x, sdlog = 0)"), + "`sdlog` must be positive", + fixed = TRUE + ) +}) diff --git a/tests/testthat/test-setup_default_parameters.R b/tests/testthat/test-setup_default_parameters.R index 5d6b6fa57..3c59d7ba9 100644 --- a/tests/testthat/test-setup_default_parameters.R +++ b/tests/testthat/test-setup_default_parameters.R @@ -23,7 +23,7 @@ test_that("`setup_default_parameters()` works with correct inputs", { colnames(result), c( "module_name", "fleet", "module_type", "label", - "age", "length", "time", "value", "estimation_type", + "age", "length", "timing", "value", "estimation_type", "distribution_type", "distribution" ) ) diff --git a/tests/testthat/test-setup_default_parameters_template.R b/tests/testthat/test-setup_default_parameters_template.R index f999875f0..c80c04cba 100644 --- a/tests/testthat/test-setup_default_parameters_template.R +++ b/tests/testthat/test-setup_default_parameters_template.R @@ -22,7 +22,7 @@ test_that("`setup_default_parameters_template()` works with correct inputs", { #' @description Test that output contains expected columns. expect_true(all(c( "module_name", "fleet", "module_type", "label", - "age", "length", "time", "value", "estimation_type", + "age", "length", "timing", "value", "estimation_type", "distribution_type", "distribution" ) %in% names(result))) }) @@ -35,7 +35,7 @@ test_that("`setup_default_parameters_template()` works with edge cases", { expect_equal(nrow(result), 0) expect_true(all(c( "module_name", "fleet", "module_type", "label", - "age", "length", "time", "value", "estimation_type", + "age", "length", "timing", "value", "estimation_type", "distribution_type", "distribution" ) %in% names(result))) }) diff --git a/tests/testthat/test-slow-parallel-caa-mle-wrappers.R b/tests/testthat/test-slow-parallel-caa-mle-wrappers.R index 148864fe0..49675428e 100644 --- a/tests/testthat/test-slow-parallel-caa-mle-wrappers.R +++ b/tests/testthat/test-slow-parallel-caa-mle-wrappers.R @@ -53,10 +53,10 @@ modified_parameters <- purrr::map(1:sim_num, \(iter_id) { tibble::tibble( fleet = "fleet1", label = "log_Fmort", - time = 1:get_n_years(data_age_length_comp), + timing = 1:get_n_years(data_age_length_comp), value = log(om_output_list[[iter_id]][["f"]]), ), - by = c("fleet", "label", "time") + by = c("fleet", "label", "timing") ) |> # Update selectivity parameters and log_q for survey1 dplyr::rows_update( @@ -67,14 +67,14 @@ modified_parameters <- purrr::map(1:sim_num, \(iter_id) { ), by = c("fleet", "label") ) |> - # Update log_devs in the Recruitment module (time steps 2-30) + # Update log_devs in the Recruitment module (timing steps 2-30) dplyr::rows_update( tibble::tibble( label = "log_devs", - time = 2:get_n_years(data_age_length_comp), + timing = 2:get_n_years(data_age_length_comp), value = om_input_list[[iter_id]][["logR.resid"]][-1] ), - by = c("label", "time") + by = c("label", "timing") ) |> # Update log_sd for log_devs in the Recruitment module dplyr::rows_update( diff --git a/vignettes/RD-re-sparsity.Rmd b/vignettes/RD-re-sparsity.Rmd index 1e2257ab7..810cfd773 100644 --- a/vignettes/RD-re-sparsity.Rmd +++ b/vignettes/RD-re-sparsity.Rmd @@ -20,7 +20,7 @@ Following statistical best practices, FIMS incorporates the use of random effect To quantify the magnitude of this penalty we investigate the performance and optimization time of a simple random effects model under equivalent sparse vs dense parameterizations. For this example we utilize two simple proxy models; an auto regressive AR1 function and a simple stock assessment model based on [babySAM](https://github.com/fishfollower/SAM/tree/master/stockassessment). These proxies should reflect the relative performance gains expected in the real-world correlated time-varying processes modeled in stock assessments such as recruitment. Each of these time-series models can be parameterized as a dense matrix, representing the traditional parameter deviations as a random-effects approach, or as a sparse matrix, as would be achieved by directly estimating annual parameter values as random effects. -These two approaches are anticipated to produce mathematically equivalent results while differing in their computational overhead. Confirming this equivalency and quantifying the runtime improvement will provide support for the change in approach relative to previous assessment models. +These two approaches are anticipated to produce mathematically equivalent results while differing in their computational overhead. Confirming this equivalency and quantifying the run-time improvement will provide support for the change in approach relative to previous assessment models. ## Methods diff --git a/vignettes/fims-demo-projections.Rmd b/vignettes/fims-demo-projections.Rmd index 43cf83bae..cc3efd1ab 100644 --- a/vignettes/fims-demo-projections.Rmd +++ b/vignettes/fims-demo-projections.Rmd @@ -100,17 +100,22 @@ data_4_projections <- data_big_with_extra_year |> dplyr::mutate( uncertainty = ifelse( type %in% c("landings") & value == -999, - 0.00999975, + "~ dlnorm(meanlog = log_landings_expected, sdlog = 0.00999975)", uncertainty ), uncertainty = ifelse( type %in% c("index") & value == -999, - 0.19804220, + "~ dlnorm(meanlog = log_index_expected, sdlog = 0.19804220)", uncertainty ), uncertainty = ifelse( - type %in% c("age_comp", "length_comp") & value == -999, - 0, + type == "age_comp" & value == -999, + "~ dmultinom(prob = agecomp_proportion, size = 0)", + uncertainty + ), + uncertainty = ifelse( + type == "length_comp" & value == -999, + "~ dmultinom(prob = lengthcomp_proportion, size = 0)", uncertainty ) ) |> @@ -136,7 +141,7 @@ parameters_projection <- setup_default_parameters(data = data_4_projections) |> tibble::tibble( fleet = "fleet1", label = "log_Fmort", - time = get_start_year(data_4_projections): + timing = get_start_year(data_4_projections): get_end_year(data_4_projections), value = log(c( 0.009459165, 0.027288858, 0.045063639, @@ -152,17 +157,17 @@ parameters_projection <- setup_default_parameters(data = data_4_projections) |> rep(0.499675368, years_of_projection) )) ), - by = c("fleet", "label", "time") + by = c("fleet", "label", "timing") ) |> # Fix the projection period log_Fmort to constant dplyr::rows_update( tibble::tibble( label = "log_Fmort", - time = (get_end_year(data_4_projections) - years_of_projection + 1): + timing = (get_end_year(data_4_projections) - years_of_projection + 1): get_end_year(data_4_projections), estimation_type = rep("constant", years_of_projection) ), - by = c("label", "time") + by = c("label", "timing") ) |> # Update selectivity parameters and log_q for survey1 dplyr::rows_update( @@ -173,11 +178,11 @@ parameters_projection <- setup_default_parameters(data = data_4_projections) |> ), by = c("fleet", "label") ) |> - # Update log_devs in the Recruitment module (time steps 2-end) + # Update log_devs in the Recruitment module (timing steps 2-end) dplyr::rows_update( tibble::tibble( label = "log_devs", - time = (get_start_year(data_4_projections) + 1): + timing = (get_start_year(data_4_projections) + 1): get_end_year(data_4_projections), value = c( 0.43787763, -0.13299042, -0.43251973, 0.64861200, 0.50640852, @@ -191,19 +196,19 @@ parameters_projection <- setup_default_parameters(data = data_4_projections) |> ), estimation_type = "random_effects" ), - by = c("label", "time") + by = c("label", "timing") ) |> # Fix the projection log recruitment deviations at zero dplyr::rows_update( tibble::tibble( label = "log_devs", - time = (get_end_year(data_4_projections) - years_of_projection + 1): + timing = (get_end_year(data_4_projections) - years_of_projection + 1): get_end_year(data_4_projections), estimation_type = rep("constant", years_of_projection), distribution_type = rep(NA, years_of_projection), distribution = rep(NA, years_of_projection) ), - by = c("label", "time") + by = c("label", "timing") ) |> # Update log_sd for log_devs in the Recruitment module dplyr::rows_update( diff --git a/vignettes/fims-demo.Rmd b/vignettes/fims-demo.Rmd index c17758f8d..849534dd6 100644 --- a/vignettes/fims-demo.Rmd +++ b/vignettes/fims-demo.Rmd @@ -141,7 +141,7 @@ parameters_4_model <- default_parameters |> tibble::tibble( fleet = "fleet1", label = "log_Fmort", - time = seq(get_n_years(data_4_model)), + timing = seq(get_n_years(data_4_model)), value = log(c( 0.009459165, 0.027288858, 0.045063639, 0.061017825, 0.048600752, 0.087420554, @@ -155,7 +155,7 @@ parameters_4_model <- default_parameters |> 0.431745298, 0.328030899, 0.499675368 )) ), - by = c("fleet", "label", "time") + by = c("fleet", "label", "timing") ) |> # Update selectivity parameters and log_q for survey1 dplyr::rows_update( @@ -166,11 +166,11 @@ parameters_4_model <- default_parameters |> ), by = c("fleet", "label") ) |> - # Update log_devs in the Recruitment module (time steps 2-30) + # Update log_devs in the Recruitment module (timing steps 2-30) dplyr::rows_update( tibble::tibble( label = "log_devs", - time = 2:get_n_years(data_4_model), + timing = 2:get_n_years(data_4_model), value = c( 0.43787763, -0.13299042, -0.43251973, 0.64861200, 0.50640852, -0.06958319, 0.30246260, -0.08257384, 0.20740372, 0.15289604, @@ -180,7 +180,7 @@ parameters_4_model <- default_parameters |> -0.19556523, 0.20094360, 0.37248740, -0.07163145 ) ), - by = c("label", "time") + by = c("label", "timing") ) |> # Update log_sd for log_devs in the Recruitment module dplyr::rows_update( diff --git a/vignettes/fims-input-data.Rmd b/vignettes/fims-input-data.Rmd index eabe117b0..3328f883e 100644 --- a/vignettes/fims-input-data.Rmd +++ b/vignettes/fims-input-data.Rmd @@ -31,7 +31,7 @@ Note: this vignette does not include anything about model configuration, which i ## Data format -The input data for a FIMS model is a single long data frame with the following columns: `r glue::glue_collapse(colnames(FIMS::data_big), sep = ", ", last = ", and ")`. Both the single table and long format are in contrast to the input format for some legacy stock assessment models where the data are provided in multiple tables, each with a different format. The SS3 data file, for instance, has separate tables for catch, indices, age compositions, length compositions, with wide tables for the composition data (each age or length bin has a different column). The use of a single, long table for data input used by FIMS makes it easier to summarize and filter the data across data types (e.g., by fleet, time step, age, etc.), as well as being better suited to the widely-used tidyverse collection of packages in R. +The input data for a FIMS model is a single long data frame with the following columns: `r glue::glue_collapse(colnames(FIMS::data_big), sep = ", ", last = ", and ")`. Both the single table and long format are in contrast to the input format for some legacy stock assessment models where the data are provided in multiple tables, each with a different format. The SS3 data file, for instance, has separate tables for catch, indices, age compositions, length compositions, with wide tables for the composition data (each age or length bin has a different column). The use of a single, long table for data input used by FIMS makes it easier to summarize and filter the data across data types (e.g., by fleet, timing, age, etc.), as well as being better suited to the widely-used tidyverse collection of packages in R. Two drawbacks of the long format are (1) some information is duplicated, and (2) the long format is harder to read and understand by looking at the raw data frame. To meet this second challenge, the FIMS team is working on functions to summarize, visualize, and check for errors in the input data. @@ -96,14 +96,28 @@ FIMS::data_big |> # time series plot of index data FIMS::data_big |> dplyr::filter(type == "index") |> - ggplot(aes(x = timing, y = value)) + - geom_point() + - geom_pointrange(aes( - ymin = qlnorm(p = 0.025, meanlog = log(value), sdlog = uncertainty), - ymax = qlnorm(p = 0.975, meanlog = log(value), sdlog = uncertainty) + ggplot2::ggplot(ggplot2::aes(x = timing, y = value)) + + ggplot2::geom_point() + + ggplot2::geom_pointrange(ggplot2::aes( + ymin = stats::qlnorm( + p = 0.025, + meanlog = log(value), + sdlog = purrr::map_dbl( + purrr::map(uncertainty, FIMS:::parse_data_distribution), + list("sdlog") + ) + ), + ymax = stats::qlnorm( + p = 0.975, + meanlog = log(value), + sdlog = purrr::map_dbl( + purrr::map(uncertainty, FIMS:::parse_data_distribution), + list("sdlog") + ) + ) )) + - geom_hline(yintercept = 0) + - ggtitle("Index time series in data_big example data set") + ggplot2::geom_hline(yintercept = 0) + + ggplot2::ggtitle("Index time series in data_big example data set") ``` ## Age compositions (`type == "age_comp"`)