Skip to contents

FIMS

The NOAA Fisheries Integrated Modeling System (FIMS) is a new modeling framework for fisheries modeling. The framework is designed to support next-generation fisheries stock assessment, ecosystem, and socioeconomic modeling. It is important to note that FIMS itself is not a model but rather a framework for creating models. The framework is made up of many modules that come together to create a model that best suits the needs of the end-user. The remainder of this vignette walks through what is absolutely necessary to run a FIMS catch-at-age model using the default settings.

Memory

Calling library(FIMS) loads the R package and Rcpp functions and modules into the R environment. The C++ code is compiled upon installation rather than loading so the call to library() should be pretty fast. Users should always run clear() prior to modeling to ensure that the C++ memory from any previous FIMS model run is cleared out.

library(FIMS)
library(ggplot2)

# clear memory
clear()

Data

Data for a FIMS model must be stored in a single data frame with a long format, e.g., data("data1", package = "FIMS"). The design is similar to running a linear model where you pass a single data frame to lm(). The long format does lead to some information being duplicated. For example, the units are listed for every row rather than stored in a single location for each data type. But, the long format facilitates using tidy functions to manipulate the data. And, a single function, i.e., FIMSFrame(), is all that is needed to prepare the data to be used in a FIMS model.

data1

A sample data frame for a catch-at-age model with both ages and lengths is stored in the package as data1. This data set is based on data that was used in Li et al. for the Model Comparison Project (github site). The length data have since been added data-raw/data1.R based on an age-length conversion matrix. See R/data1.R or ?data1 for details about the package data.

FIMSFrame()

The easiest way to prepare the data for a FIMS model is to use FIMSFrame(). This function performs several validation checks and returns an object of the S4 class called FIMSFrame. There are helper functions for working with a FIMSFrame object, e.g., get_data(), get_n_years(), get_*(). Additionally, there are helper functions for pulling data out of the S4 class in the format needed for a module, i.e., a vector, but these m_*() functions are largely meant to be used internally within the package and are only exported to allow for their use by power users wishing to manually set up.

# Bring the package data into your environment
data("data1")
# Prepare the package data for being used in a FIMS model
data_4_model <- FIMSFrame(data1)

The S4 object that we named data_4_model contains many slots (i.e., named components of the object that can be accessed) but perhaps the most interesting one is the long data frame stored in the “data” slot. The tibble stored in this slot can be accessed using get_data().

# Use show() to see what is stored in the FIMSFrame S4 class
methods::show(data_4_model)
## tbl_df of class 'FIMSFrame'
## with the following 'types': age, age-to-length-conversion, landings, length, weight-at-age, index
## # A tibble: 6 × 10
##   type  name     age length datestart dateend value unit       uncertainty  year
##   <chr> <chr>  <int>  <dbl> <date>    <date>  <dbl> <chr>            <dbl> <dbl>
## 1 age   fleet1     1     NA 1-01-01   1-12-31 0.07  proportion         200     1
## 2 age   fleet1     2     NA 1-01-01   1-12-31 0.1   proportion         200     1
## 3 age   fleet1     3     NA 1-01-01   1-12-31 0.115 proportion         200     1
## 4 age   fleet1     4     NA 1-01-01   1-12-31 0.15  proportion         200     1
## 5 age   fleet1     5     NA 1-01-01   1-12-31 0.1   proportion         200     1
## 6 age   fleet1     6     NA 1-01-01   1-12-31 0.05  proportion         200     1
## additional slots include the following:fleets:
## [1] "fleet1"  "survey1"
## n_years:
## [1] 30
## ages:
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12
## n_ages:
## [1] 12
## lengths:
##  [1]    0   50  100  150  200  250  300  350  400  450  500  550  600  650  700
## [16]  750  800  850  900  950 1000 1050 1100
## n_lengths:
## [1] 23
## start_year:
## [1] 1
## end_year:
## [1] 30
# Or, look at the structure using str()
# Increase max.level to see more of the structure
str(data_4_model, max.level = 1)
## Formal class 'FIMSFrame' [package "FIMS"] with 9 slots
# Use dplyr to subset the data for just the landings
get_data(data_4_model) |>
  dplyr::filter(type == "landings")
## # A tibble: 30 × 10
##    type   name    age length datestart dateend  value unit  uncertainty  year
##    <chr>  <chr> <int>  <dbl> <date>    <date>   <dbl> <chr>       <dbl> <dbl>
##  1 landi… flee…    NA     NA 1-01-01   1-12-31   162. mt         0.0100     1
##  2 landi… flee…    NA     NA 2-01-01   2-12-31   461. mt         0.0100     2
##  3 landi… flee…    NA     NA 3-01-01   3-12-31   747. mt         0.0100     3
##  4 landi… flee…    NA     NA 4-01-01   4-12-31   997. mt         0.0100     4
##  5 landi… flee…    NA     NA 5-01-01   5-12-31   768. mt         0.0100     5
##  6 landi… flee…    NA     NA 6-01-01   6-12-31  1344. mt         0.0100     6
##  7 landi… flee…    NA     NA 7-01-01   7-12-31  1319. mt         0.0100     7
##  8 landi… flee…    NA     NA 8-01-01   8-12-31  2598. mt         0.0100     8
##  9 landi… flee…    NA     NA 9-01-01   9-12-31  1426. mt         0.0100     9
## 10 landi… flee…    NA     NA 10-01-01  10-12-31 1644. mt         0.0100    10
## # ℹ 20 more rows

The data contains the following fleets:

  • A single fishery fleet with age- and length-composition, weight-at-age, and landings data
  • A single survey with age- and length-composition and index data

Parameters

The parameters that are in the model will depend on which modules are used from the FIMS framework. This combination of modules rather than the use of a control file negates the need for complicated if{} else{} statements in the code.

create_default_parameters()

By passing the data to create_default_parameters() the function can tailor the defaults based on how many fleets there are and what data types exist. For example, if you have three fleets, then create_default_parameters() will set up three logistic selectivity modules.

Modules that are available in FIMS are known as reference classes in the C++ code. Each reference class acts as an interface between R and the underlining C++ code that defines FIMS. Several reference classes exist and several more will be created in the future. The beauty of having modules rather than a control file really comes out when more reference classes are created because each reference class can be accessed through R by itself to build up a model rather than needing to modify a control file for future features.

By just passing lists of fleet specifications and the data to create_default_parameters(), the default values for parameters that relate to fleet(s), recruitment, growth, and maturity modules can be created. For example,

  • “BevertonHoltRecruitment” for the recruitment module
  • “DnormDistribution” for recruitment deviations (log_devs)
  • “EWAAgrowth” for the growth module, and
  • “LogisticMaturity” for maturity module.
# Define the same fleet specifications for fleet1 and survey1
fleet1 <- survey1 <- list(
  selectivity = list(form = "LogisticSelectivity"),
  data_distribution = c(
    Landings = "DlnormDistribution",
    Index = "DlnormDistribution",
    AgeComp = "DmultinomDistribution",
    LengthComp = "DmultinomDistribution"
  )
)

# Create default recruitment, growth, and maturity parameters
default_parameters <- data_4_model |>
  create_default_parameters(
    fleets = list(fleet1 = fleet1, survey1 = survey1)
  )

The argument names and their corresponding default values of the create_default_parameters() function can be displayed using args().

args(create_default_parameters)
## function (data, fleets, recruitment = list(form = "BevertonHoltRecruitment", 
##     process_distribution = c(log_devs = "DnormDistribution")), 
##     growth = list(form = "EWAAgrowth"), maturity = list(form = "LogisticMaturity")) 
## NULL

update_parameters()

In the future, the developers of FIMS may update the default parameters to experiment with different values. Regardless, you can still use create_default_parameters() as a starting point because it will provide information on the appropriate dimensions and necessary elements that the final list must contain. And, although it is a good idea to modify the returned defaults with update_parameters(), you can update the list manually.

In the code below, update_parameters() is used to adjust the fishing mortality, selectivity, maturity, and population parameters from their default values. If the parameters are estimated, these updates will change their starting values, and if they are fixed, these updates will change their value used in the model.

# Each call to update_parameters() returns the full list so the pipe can be
# used to daisy chain all of these updates together to a new object called
# parameters that will be used to fit the model
parameters <- default_parameters |>
  update_parameters(
    modified_parameters = list(
      fleet1 = list(
        Fleet.log_Fmort.value = log(c(
          0.009459165, 0.027288858, 0.045063639,
          0.061017825, 0.048600752, 0.087420554,
          0.088447204, 0.186607929, 0.109008958,
          0.132704335, 0.150615473, 0.161242955,
          0.116640187, 0.169346119, 0.180191913,
          0.161240483, 0.314573212, 0.257247574,
          0.254887252, 0.251462108, 0.349101406,
          0.254107720, 0.418478117, 0.345721184,
          0.343685540, 0.314171227, 0.308026829,
          0.431745298, 0.328030899, 0.499675368
        ))
      )
    )
  ) |>
  update_parameters(
    modified_parameters = list(
      survey1 = list(
        LogisticSelectivity.inflection_point.value = 1.5,
        LogisticSelectivity.slope.value = 2,
        Fleet.log_q.value = log(3.315143e-07)
      )
    )
  ) |>
  update_parameters(
    modified_parameters = list(
      recruitment = list(
        BevertonHoltRecruitment.log_rzero.value = log(1e+06),
        BevertonHoltRecruitment.log_devs.value = c(
          0.43787763, -0.13299042, -0.43251973, 0.64861200, 0.50640852,
          -0.06958319, 0.30246260, -0.08257384, 0.20740372, 0.15289604,
          -0.21709207, -0.13320626, 0.11225374, -0.10650836, 0.26877132,
          0.24094126, -0.54480751, -0.23680557, -0.58483386, 0.30122785,
          0.21930545, -0.22281699, -0.51358369, 0.15740234, -0.53988240,
          -0.19556523, 0.20094360, 0.37248740, -0.07163145
        ),
        BevertonHoltRecruitment.log_devs.estimation_type = "fixed_effects"
      )
    )
  ) |>
  update_parameters(
    modified_parameters = list(
      maturity = list(
        LogisticMaturity.inflection_point.value = 2.25,
        LogisticMaturity.inflection_point.estimation_type = "constant",
        LogisticMaturity.slope.value = 3,
        LogisticMaturity.slope.estimation_type = "constant"
      )
    )
  ) |>
  update_parameters(
    modified_parameters = list(
      population = list(
        Population.log_init_naa.value = c(
          13.80944, 13.60690, 13.40217, 13.19525, 12.98692, 12.77791,
          12.56862, 12.35922, 12.14979, 11.94034, 11.73088, 13.18755
        )
      )
    )
  )

purrr::map_vec() can be used to compare the length of individual updated parameter vectors from parameters with the length of individual parameter vectors from default_parameters.

default_fleet1 <- purrr::map_vec(
  default_parameters[["parameters"]][["fleet1"]],
  \(x) length(x)
)
updated_fleet1 <- purrr::map_vec(
  parameters[["parameters"]][["fleet1"]],
  \(x) length(x)
)
knitr::kable(cbind(default_fleet1, updated_fleet1))
default_fleet1 updated_fleet1
LogisticSelectivity.inflection_point.value 1 1
LogisticSelectivity.inflection_point.estimation_type 1 1
LogisticSelectivity.slope.value 1 1
LogisticSelectivity.slope.estimation_type 1 1
Fleet.log_q.value 1 1
Fleet.log_q.estimation_type 1 1
Fleet.log_Fmort.value 30 30
Fleet.log_Fmort.estimation_type 1 1
DlnormDistribution.log_sd.value 30 30
DlnormDistribution.log_sd.estimation_type 1 1

purrr::map_df() can be used to summarize parameter vector lengths across all modules.

purrr::map_df(
  parameters[["parameters"]], \(x) purrr::map_vec(x, length),
  .id = "module"
) |>
  tibble::column_to_rownames(var = "module") |>
  t()
##                                                                  fleet1 survey1
## LogisticSelectivity.inflection_point.value                            1       1
## LogisticSelectivity.inflection_point.estimation_type                  1       1
## LogisticSelectivity.slope.value                                       1       1
## LogisticSelectivity.slope.estimation_type                             1       1
## Fleet.log_q.value                                                     1       1
## Fleet.log_q.estimation_type                                           1       1
## Fleet.log_Fmort.value                                                30      30
## Fleet.log_Fmort.estimation_type                                       1       1
## DlnormDistribution.log_sd.value                                      30      30
## DlnormDistribution.log_sd.estimation_type                             1       1
## BevertonHoltRecruitment.log_rzero.value                              NA      NA
## BevertonHoltRecruitment.log_rzero.estimation_type                    NA      NA
## BevertonHoltRecruitment.log_r.value                                  NA      NA
## BevertonHoltRecruitment.log_r.estimation_type                        NA      NA
## BevertonHoltRecruitment.logit_steep.value                            NA      NA
## BevertonHoltRecruitment.logit_steep.estimation_type                  NA      NA
## BevertonHoltRecruitment.log_devs.value                               NA      NA
## BevertonHoltRecruitment.log_devs.estimation_type                     NA      NA
## BevertonHoltRecruitment.log_expected_recruitment.value               NA      NA
## BevertonHoltRecruitment.log_expected_recruitment.estimation_type     NA      NA
## DnormDistribution.log_sd.value                                       NA      NA
## DnormDistribution.log_sd.estimation_type                             NA      NA
## DnormDistribution.x.value                                            NA      NA
## DnormDistribution.x.estimation_type                                  NA      NA
## DnormDistribution.expected_values.value                              NA      NA
## DnormDistribution.expected_values.estimation_type                    NA      NA
## LogisticMaturity.inflection_point.value                              NA      NA
## LogisticMaturity.inflection_point.estimation_type                    NA      NA
## LogisticMaturity.slope.value                                         NA      NA
## LogisticMaturity.slope.estimation_type                               NA      NA
## Population.log_M.value                                               NA      NA
## Population.log_M.estimation_type                                     NA      NA
## Population.log_init_naa.value                                        NA      NA
## Population.log_init_naa.estimation_type                              NA      NA
##                                                                  recruitment
## LogisticSelectivity.inflection_point.value                                NA
## LogisticSelectivity.inflection_point.estimation_type                      NA
## LogisticSelectivity.slope.value                                           NA
## LogisticSelectivity.slope.estimation_type                                 NA
## Fleet.log_q.value                                                         NA
## Fleet.log_q.estimation_type                                               NA
## Fleet.log_Fmort.value                                                     NA
## Fleet.log_Fmort.estimation_type                                           NA
## DlnormDistribution.log_sd.value                                           NA
## DlnormDistribution.log_sd.estimation_type                                 NA
## BevertonHoltRecruitment.log_rzero.value                                    1
## BevertonHoltRecruitment.log_rzero.estimation_type                          1
## BevertonHoltRecruitment.log_r.value                                       29
## BevertonHoltRecruitment.log_r.estimation_type                              1
## BevertonHoltRecruitment.logit_steep.value                                  1
## BevertonHoltRecruitment.logit_steep.estimation_type                        1
## BevertonHoltRecruitment.log_devs.value                                    29
## BevertonHoltRecruitment.log_devs.estimation_type                           1
## BevertonHoltRecruitment.log_expected_recruitment.value                    31
## BevertonHoltRecruitment.log_expected_recruitment.estimation_type           1
## DnormDistribution.log_sd.value                                             1
## DnormDistribution.log_sd.estimation_type                                   1
## DnormDistribution.x.value                                                 30
## DnormDistribution.x.estimation_type                                        1
## DnormDistribution.expected_values.value                                   30
## DnormDistribution.expected_values.estimation_type                          1
## LogisticMaturity.inflection_point.value                                   NA
## LogisticMaturity.inflection_point.estimation_type                         NA
## LogisticMaturity.slope.value                                              NA
## LogisticMaturity.slope.estimation_type                                    NA
## Population.log_M.value                                                    NA
## Population.log_M.estimation_type                                          NA
## Population.log_init_naa.value                                             NA
## Population.log_init_naa.estimation_type                                   NA
##                                                                  maturity
## LogisticSelectivity.inflection_point.value                             NA
## LogisticSelectivity.inflection_point.estimation_type                   NA
## LogisticSelectivity.slope.value                                        NA
## LogisticSelectivity.slope.estimation_type                              NA
## Fleet.log_q.value                                                      NA
## Fleet.log_q.estimation_type                                            NA
## Fleet.log_Fmort.value                                                  NA
## Fleet.log_Fmort.estimation_type                                        NA
## DlnormDistribution.log_sd.value                                        NA
## DlnormDistribution.log_sd.estimation_type                              NA
## BevertonHoltRecruitment.log_rzero.value                                NA
## BevertonHoltRecruitment.log_rzero.estimation_type                      NA
## BevertonHoltRecruitment.log_r.value                                    NA
## BevertonHoltRecruitment.log_r.estimation_type                          NA
## BevertonHoltRecruitment.logit_steep.value                              NA
## BevertonHoltRecruitment.logit_steep.estimation_type                    NA
## BevertonHoltRecruitment.log_devs.value                                 NA
## BevertonHoltRecruitment.log_devs.estimation_type                       NA
## BevertonHoltRecruitment.log_expected_recruitment.value                 NA
## BevertonHoltRecruitment.log_expected_recruitment.estimation_type       NA
## DnormDistribution.log_sd.value                                         NA
## DnormDistribution.log_sd.estimation_type                               NA
## DnormDistribution.x.value                                              NA
## DnormDistribution.x.estimation_type                                    NA
## DnormDistribution.expected_values.value                                NA
## DnormDistribution.expected_values.estimation_type                      NA
## LogisticMaturity.inflection_point.value                                 1
## LogisticMaturity.inflection_point.estimation_type                       1
## LogisticMaturity.slope.value                                            1
## LogisticMaturity.slope.estimation_type                                  1
## Population.log_M.value                                                 NA
## Population.log_M.estimation_type                                       NA
## Population.log_init_naa.value                                          NA
## Population.log_init_naa.estimation_type                                NA
##                                                                  population
## LogisticSelectivity.inflection_point.value                               NA
## LogisticSelectivity.inflection_point.estimation_type                     NA
## LogisticSelectivity.slope.value                                          NA
## LogisticSelectivity.slope.estimation_type                                NA
## Fleet.log_q.value                                                        NA
## Fleet.log_q.estimation_type                                              NA
## Fleet.log_Fmort.value                                                    NA
## Fleet.log_Fmort.estimation_type                                          NA
## DlnormDistribution.log_sd.value                                          NA
## DlnormDistribution.log_sd.estimation_type                                NA
## BevertonHoltRecruitment.log_rzero.value                                  NA
## BevertonHoltRecruitment.log_rzero.estimation_type                        NA
## BevertonHoltRecruitment.log_r.value                                      NA
## BevertonHoltRecruitment.log_r.estimation_type                            NA
## BevertonHoltRecruitment.logit_steep.value                                NA
## BevertonHoltRecruitment.logit_steep.estimation_type                      NA
## BevertonHoltRecruitment.log_devs.value                                   NA
## BevertonHoltRecruitment.log_devs.estimation_type                         NA
## BevertonHoltRecruitment.log_expected_recruitment.value                   NA
## BevertonHoltRecruitment.log_expected_recruitment.estimation_type         NA
## DnormDistribution.log_sd.value                                           NA
## DnormDistribution.log_sd.estimation_type                                 NA
## DnormDistribution.x.value                                                NA
## DnormDistribution.x.estimation_type                                      NA
## DnormDistribution.expected_values.value                                  NA
## DnormDistribution.expected_values.estimation_type                        NA
## LogisticMaturity.inflection_point.value                                  NA
## LogisticMaturity.inflection_point.estimation_type                        NA
## LogisticMaturity.slope.value                                             NA
## LogisticMaturity.slope.estimation_type                                   NA
## Population.log_M.value                                                  360
## Population.log_M.estimation_type                                          1
## Population.log_init_naa.value                                            12
## Population.log_init_naa.estimation_type                                   1

Fit

With data and parameters in place, we can now initialize modules using initialize_fims() and fit the model using fit_fims().

initialize_fims()

The list returned by create_default_parameters() has two elements, parameters and modules. But, these are just lists of lists containing specifications. Nothing has been created in memory as of yet. To actually initialize the modules specified in parameters[["modules"]], initialize_fims() needs to be called. This function takes all of the specifications and matches them with the appropriate data to initialize a module and create the pointers to the memory.

fit_fims()

The list returned from initialize_fims() can be passed to the parameter of fit_fims() called input to run a FIMS model. If optimize = FALSE, the model will not actually be optimized but instead just checked to ensure it is a viable model. When optimize = TRUE, the model will be fit using stats::nlminb() and an object of the class FIMSFit will be returned.

Example

# Run the model without optimization to help ensure a viable model
test_fit <- parameters |>
  initialize_fims(data = data_4_model) |>
  fit_fims(optimize = FALSE)
clear()
# Run the  model with optimization
fit <- parameters |>
  initialize_fims(data = data_4_model) |>
  fit_fims(optimize = TRUE)
##  Starting optimization ...
##  Restarting optimizer 3 times to improve gradient.
##  Maximum gradient went from 0.00739 to 0.00331 after 3 steps.
##  Finished optimization
##  Finished sdreport
##  FIMS model version: 0.5.0
##  Total run time was 6.51132 seconds
##  Number of parameters: total=77, fixed_effects=77, and random_effects=0
##  Maximum gradient= 0.00331
##  Negative log likelihood (NLL):
##  Marginal NLL= 3239.98416
##  Total NLL= 3239.98416
##  Terminal SB= 1853.25812
clear()
# finalize is not working, you get the following error when you do this code
# Error: basic_string::_M_create
# if you run clear() ahead of the following line the error will go away
# fit_output <- finalize(fit@opt$par, fit@obj$fn, fit@obj$gr)

# Get information about the model and print a few characters to the screen
recruitment_log <- get_log_module("information")
substr(recruitment_log, 1, 100)
## [1] "[\n]"

The results can be plotted with either base R or using {ggplot2}.

index_results <- data.frame(
  observed = m_index(data_4_model, "survey1"),
  expected = get_report(fit)[["index_exp"]][[2]],
  years = get_start_year(data_4_model):get_end_year(data_4_model)
)
print(index_results)
##       observed    expected years
## 1  0.006117418 0.006522883     1
## 2  0.007156588 0.006583769     2
## 3  0.006553376 0.006588869     3
## 4  0.006501745 0.006469924     4
## 5  0.006562572 0.006230236     5
## 6  0.006467796 0.006199557     6
## 7  0.005248220 0.006046247     7
## 8  0.006568311 0.005874548     8
## 9  0.005516030 0.005245321     9
## 10 0.005614254 0.005042116    10
## 11 0.004269091 0.004761138    11
## 12 0.004842784 0.004432323    12
## 13 0.004018395 0.004078090    13
## 14 0.004741539 0.003928076    14
## 15 0.003444776 0.003655184    15
## 16 0.004062555 0.003424690    16
## 17 0.003790370 0.003322162    17
## 18 0.002892876 0.002873339    18
## 19 0.002765623 0.002600646    19
## 20 0.002349581 0.002335867    20
## 21 0.002007219 0.002125365    21
## 22 0.001883646 0.001897374    22
## 23 0.001908972 0.001943536    23
## 24 0.001596426 0.001722025    24
## 25 0.001248675 0.001585901    25
## 26 0.001459593 0.001490177    26
## 27 0.001255182 0.001445525    27
## 28 0.001314476 0.001423363    28
## 29 0.001364266 0.001353681    29
## 30 0.001343069 0.001455624    30
ggplot2::ggplot(index_results, ggplot2::aes(x = years, y = observed)) +
  ggplot2::geom_point() +
  ggplot2::xlab("Year") +
  ggplot2::ylab("Index (mt)") +
  ggplot2::geom_line(ggplot2::aes(x = years, y = expected), color = "blue") +
  ggplot2::theme_bw()

landings_results <- data.frame(
  observed = m_landings(data_4_model, fleet = "fleet1"),
  expected = get_report(fit)[["landings_exp"]][[1]],
  years = get_start_year(data_4_model):get_end_year(data_4_model)
)
print(landings_results)
##     observed  expected years
## 1   161.6455  161.6447     1
## 2   461.0895  461.0651     2
## 3   747.2900  747.2166     3
## 4   996.9711  996.9583     4
## 5   767.5477  767.5261     5
## 6  1343.8609 1343.6351     6
## 7  1319.0663 1318.7179     7
## 8  2597.5051 2595.9616     8
## 9  1425.8608 1425.2560     9
## 10 1643.6742 1642.6459    10
## 11 1770.8245 1769.3562    11
## 12 1751.2628 1749.9570    12
## 13 1194.0655 1193.5915    13
## 14 1638.4036 1637.7572    14
## 15 1584.1855 1584.0130    15
## 16 1333.1820 1332.9450    16
## 17 2325.4425 2324.3325    17
## 18 1693.7841 1693.3032    18
## 19 1531.8238 1531.8923    19
## 20 1358.8162 1360.2285    20
## 21 1638.8798 1640.7833    21
## 22 1077.3713 1077.5113    22
## 23 1647.9104 1647.3195    23
## 24 1224.4634 1225.2509    24
## 25 1126.9182 1127.2119    25
## 26  956.4529  956.9399    26
## 27  902.1621  903.0091    27
## 28 1168.8217 1170.4514    28
## 29  852.3819  852.3091    29
## 30 1274.3777 1274.3969    30
ggplot2::ggplot(landings_results, ggplot2::aes(x = years, y = observed)) +
  ggplot2::geom_point() +
  ggplot2::xlab("Year") +
  ggplot2::ylab("Index (mt)") +
  ggplot2::geom_line(ggplot2::aes(x = years, y = expected), color = "blue") +
  ggplot2::theme_bw()

Sensitivities

Multiple fits, i.e., sensitivity runs, can be set up by modifying the parameter list using update_parameters() or changing the data that is used to fit the model.

Initial values

For example, one could change the initial value used for the slope of the logistic curve for the survey to see if the terminal estimate changes due to changes to the initial value.

parameters_high_slope <- parameters |>
  update_parameters(
    modified_parameters = list(
      survey1 = list(
        LogisticSelectivity.slope.value = 2.5
      )
    )
  )

parameters_low_slope <- parameters |>
  update_parameters(
    modified_parameters = list(
      survey1 = list(
        LogisticSelectivity.slope.value = 1
      )
    )
  )

high_slope_fit <- parameters_high_slope |>
  initialize_fims(data = data_4_model) |>
  fit_fims(optimize = TRUE)
##  Starting optimization ...
##  Restarting optimizer 3 times to improve gradient.
##  Maximum gradient went from 0.07023 to 0.00122 after 3 steps.
##  Finished optimization
##  Finished sdreport
##  FIMS model version: 0.5.0
##  Total run time was 6.36162 seconds
##  Number of parameters: total=77, fixed_effects=77, and random_effects=0
##  Maximum gradient= 0.00122
##  Negative log likelihood (NLL):
##  Marginal NLL= 3239.98416
##  Total NLL= 3239.98416
##  Terminal SB= 1853.26753
clear()

low_slope_fit <- parameters_low_slope |>
  initialize_fims(data = data_4_model) |>
  fit_fims(optimize = TRUE)
##  Starting optimization ...
##  Restarting optimizer 3 times to improve gradient.
##  Maximum gradient went from 0.01168 to 0.00157 after 3 steps.
##  Finished optimization
##  Finished sdreport
##  FIMS model version: 0.5.0
##  Total run time was 6.35972 seconds
##  Number of parameters: total=77, fixed_effects=77, and random_effects=0
##  Maximum gradient= 0.00157
##  Negative log likelihood (NLL):
##  Marginal NLL= 3239.98416
##  Total NLL= 3239.98416
##  Terminal SB= 1853.25978
clear()

Age only

The same model can be fit to just the age data, removing the length data.

# Define fleet and survey with age-specific data distribution
fleet1 <- survey1 <- list(
  selectivity = list(form = "LogisticSelectivity"),
  data_distribution = c(
    Index = "DlnormDistribution",
    Landings = "DlnormDistribution",
    AgeComp = "DmultinomDistribution"
  )
)

# Create default parameters, update with modified values, initialize FIMS,
# and fit the model
age_only_fit <- data_4_model |>
  create_default_parameters(
    fleets = list(fleet1 = fleet1, survey1 = survey1)
  ) |>
  # update_parameters(modified_parameters = parameters$parameters) |>
  initialize_fims(data = data_4_model) |>
  fit_fims(optimize = TRUE)
##  Starting optimization ...
##  Restarting optimizer 3 times to improve gradient.
##  Maximum gradient went from 0.01636 to 0.001 after 3 steps.
##  Finished optimization
##  Finished sdreport
##  FIMS model version: 0.5.0
##  Total run time was 0.59721 seconds
##  Number of parameters: total=79, fixed_effects=79, and random_effects=0
##  Maximum gradient= 0.001
##  Negative log likelihood (NLL):
##  Marginal NLL= 1619.91297
##  Total NLL= 1619.91297
##  Terminal SB= 387.938
clear()

Length

The same model can be fit to just the length data, removing the age data.

# Define fleet and survey with length-specific data distribution
fleet1 <- survey1 <- list(
  selectivity = list(form = "LogisticSelectivity"),
  data_distribution = c(
    Index = "DlnormDistribution",
    Landings = "DlnormDistribution",
    LengthComp = "DmultinomDistribution"
  )
)

# Create default parameters, update with modified values, initialize FIMS,
# and fit the model
length_only_fit <- data_4_model |>
  create_default_parameters(
    fleets = list(fleet1 = fleet1, survey1 = survey1)
  ) |>
  update_parameters(modified_parameters = parameters$parameters) |>
  initialize_fims(data = data_4_model) |>
  fit_fims(optimize = TRUE)
##  Starting optimization ...
##  Restarting optimizer 3 times to improve gradient.
##  Maximum gradient went from 0.01735 to 0.00018 after 3 steps.
##  Finished optimization
##  Finished sdreport
##  FIMS model version: 0.5.0
##  Total run time was 4.89931 seconds
##  Number of parameters: total=77, fixed_effects=77, and random_effects=0
##  Maximum gradient= 0.00018
##  Negative log likelihood (NLL):
##  Marginal NLL= 1551.78735
##  Total NLL= 1551.78735
##  Terminal SB= 1778.81313
# cat(fit_output)
clear()