Introducing the Fisheries Integrated Modeling System (FIMS)
Source:vignettes/fims-demo.Rmd
fims-demo.Rmd
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.
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-to-length-conversion, age_comp, landings, length_comp, weight-at-age, index
## # A tibble: 6 × 10
## type name age length datestart dateend value unit uncertainty
## <chr> <chr> <int> <dbl> <date> <date> <dbl> <chr> <dbl>
## 1 age-to-le… flee… 1 0 1-01-01 1-12-31 1.26e-16 proporti… 200
## 2 age-to-le… flee… 1 50 1-01-01 1-12-31 8.39e-11 proporti… 200
## 3 age-to-le… flee… 1 100 1-01-01 1-12-31 2.30e- 6 proporti… 200
## 4 age-to-le… flee… 1 150 1-01-01 1-12-31 2.74e- 3 proporti… 200
## 5 age-to-le… flee… 1 200 1-01-01 1-12-31 1.63e- 1 proporti… 200
## 6 age-to-le… flee… 1 250 1-01-01 1-12-31 6.32e- 1 proporti… 200
## # ℹ 1 more variable: year <dbl>
## 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.01000 1
## 2 landi… flee… NA NA 2-01-01 2-12-31 461. mt 0.01000 2
## 3 landi… flee… NA NA 3-01-01 3-12-31 747. mt 0.01000 3
## 4 landi… flee… NA NA 4-01-01 4-12-31 997. mt 0.01000 4
## 5 landi… flee… NA NA 5-01-01 5-12-31 768. mt 0.01000 5
## 6 landi… flee… NA NA 6-01-01 6-12-31 1344. mt 0.01000 6
## 7 landi… flee… NA NA 7-01-01 7-12-31 1319. mt 0.01000 7
## 8 landi… flee… NA NA 8-01-01 8-12-31 2598. mt 0.01000 8
## 9 landi… flee… NA NA 9-01-01 9-12-31 1426. mt 0.01000 9
## 10 landi… flee… NA NA 10-01-01 10-12-31 1644. mt 0.01000 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
configurations
create_default_configurations()
The create_default_configurations()
function is designed
to generate a set of default configurations for the various components
of a FIMS model. This includes configurations for fleets, growth,
maturity, and recruitment modules. By leveraging the structure of the
input data, the function can automatically set up initial configurations
for each module. By passing the data and configurations 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_configurations()
will set up three logistic
selectivity modules.
# Create default configurations based on the data
default_configurations <- create_default_configurations(data = data_4_model) |>
print()
##
[38;5;246m# A tibble: 7 × 4
[39m
## model_family module_name fleet_name data
##
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<list>
[39m
[23m
##
[38;5;250m1
[39m catch_at_age Data fleet1
[38;5;246m<tibble [3 × 4]>
[39m
##
[38;5;250m2
[39m catch_at_age Selectivity fleet1
[38;5;246m<tibble [1 × 4]>
[39m
##
[38;5;250m3
[39m catch_at_age Data survey1
[38;5;246m<tibble [3 × 4]>
[39m
##
[38;5;250m4
[39m catch_at_age Selectivity survey1
[38;5;246m<tibble [1 × 4]>
[39m
##
[38;5;250m5
[39m catch_at_age Growth
[31mNA
[39m
[38;5;246m<tibble [1 × 4]>
[39m
##
[38;5;250m6
[39m catch_at_age Maturity
[31mNA
[39m
[38;5;246m<tibble [1 × 4]>
[39m
##
[38;5;250m7
[39m catch_at_age Recruitment
[31mNA
[39m
[38;5;246m<tibble [1 × 4]>
[39m
# The output is a nested tibble, with details in the `data` column.
default_configurations_unnested <- default_configurations |>
tidyr::unnest(cols = data) |>
print()
##
[38;5;246m# A tibble: 11 × 7
[39m
## model_family module_name fleet_name module_type distribution_link
##
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
##
[38;5;250m 1
[39m catch_at_age Data fleet1 AgeComp AgeComp
##
[38;5;250m 2
[39m catch_at_age Data fleet1 Landings Landings
##
[38;5;250m 3
[39m catch_at_age Data fleet1 LengthComp LengthComp
##
[38;5;250m 4
[39m catch_at_age Selectivity fleet1 Logistic
[31mNA
[39m
##
[38;5;250m 5
[39m catch_at_age Data survey1 AgeComp AgeComp
##
[38;5;250m 6
[39m catch_at_age Data survey1 Index Index
##
[38;5;250m 7
[39m catch_at_age Data survey1 LengthComp LengthComp
##
[38;5;250m 8
[39m catch_at_age Selectivity survey1 Logistic
[31mNA
[39m
##
[38;5;250m 9
[39m catch_at_age Growth
[31mNA
[39m EWAA
[31mNA
[39m
##
[38;5;250m10
[39m catch_at_age Maturity
[31mNA
[39m Logistic
[31mNA
[39m
##
[38;5;250m11
[39m catch_at_age Recruitment
[31mNA
[39m BevertonHolt log_devs
##
[38;5;246m# ℹ 2 more variables: distribution_type <chr>, distribution <chr>
[39m
Update configurations
The default_configurations are just a starting point. Functions
(e.g., rows_*()
) from dplyr can be used to
modify the default configurations as needed. For example, logistic
selectivity for survey1 can be changed to double logistic
selectivity.
# Update the module_type for survey1's selectivity
updated_configurations <- default_configurations_unnested |>
dplyr::rows_update(
y = tibble::tibble(
module_name = c("Selectivity"),
fleet_name = c("survey1"),
module_type = c("DoubleLogistic")
),
by = c("module_name", "fleet_name")
) |>
print()
##
[38;5;246m# A tibble: 11 × 7
[39m
## model_family module_name fleet_name module_type distribution_link
##
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
##
[38;5;250m 1
[39m catch_at_age Data fleet1 AgeComp AgeComp
##
[38;5;250m 2
[39m catch_at_age Data fleet1 Landings Landings
##
[38;5;250m 3
[39m catch_at_age Data fleet1 LengthComp LengthComp
##
[38;5;250m 4
[39m catch_at_age Selectivity fleet1 Logistic
[31mNA
[39m
##
[38;5;250m 5
[39m catch_at_age Data survey1 AgeComp AgeComp
##
[38;5;250m 6
[39m catch_at_age Data survey1 Index Index
##
[38;5;250m 7
[39m catch_at_age Data survey1 LengthComp LengthComp
##
[38;5;250m 8
[39m catch_at_age Selectivity survey1 DoubleLogistic
[31mNA
[39m
##
[38;5;250m 9
[39m catch_at_age Growth
[31mNA
[39m EWAA
[31mNA
[39m
##
[38;5;250m10
[39m catch_at_age Maturity
[31mNA
[39m Logistic
[31mNA
[39m
##
[38;5;250m11
[39m catch_at_age Recruitment
[31mNA
[39m BevertonHolt log_devs
##
[38;5;246m# ℹ 2 more variables: distribution_type <chr>, distribution <chr>
[39m
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()
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 the configurations 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,
- “BevertonHolt” for the recruitment module
- “Dnorm” distribution for recruitment deviations (log_devs)
- “EWAA” for the Growth module, and
- “Logistic” for Maturity module.
# Create default parameters based on default_configurations and data
default_parameters <- create_default_parameters(
configurations = default_configurations,
data = data_4_model
) |>
print()
##
[38;5;246m# A tibble: 10 × 4
[39m
## model_family module_name fleet_name data
##
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<list>
[39m
[23m
##
[38;5;250m 1
[39m catch_at_age Selectivity fleet1
[38;5;246m<tibble [2 × 10]>
[39m
##
[38;5;250m 2
[39m catch_at_age Fleet fleet1
[38;5;246m<tibble [31 × 10]>
[39m
##
[38;5;250m 3
[39m catch_at_age Data fleet1
[38;5;246m<tibble [32 × 10]>
[39m
##
[38;5;250m 4
[39m catch_at_age Selectivity survey1
[38;5;246m<tibble [2 × 10]>
[39m
##
[38;5;250m 5
[39m catch_at_age Fleet survey1
[38;5;246m<tibble [31 × 10]>
[39m
##
[38;5;250m 6
[39m catch_at_age Data survey1
[38;5;246m<tibble [32 × 10]>
[39m
##
[38;5;250m 7
[39m catch_at_age Recruitment
[31mNA
[39m
[38;5;246m<tibble [150 × 10]>
[39m
##
[38;5;250m 8
[39m catch_at_age Maturity
[31mNA
[39m
[38;5;246m<tibble [2 × 10]>
[39m
##
[38;5;250m 9
[39m catch_at_age Population
[31mNA
[39m
[38;5;246m<tibble [372 × 10]>
[39m
##
[38;5;250m10
[39m catch_at_age Growth
[31mNA
[39m
[38;5;246m<tibble [1 × 10]>
[39m
# Unnest the default_parameters to see the detailed parameters
default_parameters_unnested <- tidyr::unnest(default_parameters, cols = data) |>
print()
##
[38;5;246m# A tibble: 655 × 13
[39m
## model_family module_name fleet_name module_type label distribution_link age
##
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
##
[38;5;250m 1
[39m catch_at_age Selectivity fleet1 Logistic infl…
[31mNA
[39m
[31mNA
[39m
##
[38;5;250m 2
[39m catch_at_age Selectivity fleet1 Logistic slope
[31mNA
[39m
[31mNA
[39m
##
[38;5;250m 3
[39m catch_at_age Fleet fleet1
[31mNA
[39m log_q
[31mNA
[39m
[31mNA
[39m
##
[38;5;250m 4
[39m catch_at_age Fleet fleet1
[31mNA
[39m log_…
[31mNA
[39m
[31mNA
[39m
##
[38;5;250m 5
[39m catch_at_age Fleet fleet1
[31mNA
[39m log_…
[31mNA
[39m
[31mNA
[39m
##
[38;5;250m 6
[39m catch_at_age Fleet fleet1
[31mNA
[39m log_…
[31mNA
[39m
[31mNA
[39m
##
[38;5;250m 7
[39m catch_at_age Fleet fleet1
[31mNA
[39m log_…
[31mNA
[39m
[31mNA
[39m
##
[38;5;250m 8
[39m catch_at_age Fleet fleet1
[31mNA
[39m log_…
[31mNA
[39m
[31mNA
[39m
##
[38;5;250m 9
[39m catch_at_age Fleet fleet1
[31mNA
[39m log_…
[31mNA
[39m
[31mNA
[39m
##
[38;5;250m10
[39m catch_at_age Fleet fleet1
[31mNA
[39m log_…
[31mNA
[39m
[31mNA
[39m
##
[38;5;246m# ℹ 645 more rows
[39m
##
[38;5;246m# ℹ 6 more variables: length <dbl>, time <dbl>, value <dbl>,
[39m
##
[38;5;246m# estimation_type <chr>, distribution_type <chr>, distribution <chr>
[39m
Update parameters
Functions (e.g., rows_*()
) from dplyr can
be used to update the default parameters as needed.
In the code below, rows_update()
is used to adjust the
fishing mortality, selectivity, maturity, and population parameters from
their default values.
parameters_4_model <- default_parameters |>
tidyr::unnest(cols = data) |>
# Update log_Fmort initial values for Fleet1
dplyr::rows_update(
tibble::tibble(
fleet_name = "fleet1",
label = "log_Fmort",
time = 1:30,
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
))
),
by = c("fleet_name", "label", "time")
) |>
# Update selectivity parameters and log_q for survey1
dplyr::rows_update(
tibble::tibble(
fleet_name = "survey1",
label = c("inflection_point", "slope", "log_q"),
value = c(1.5, 2, log(3.315143e-07))
),
by = c("fleet_name", "label")
) |>
# Update log_devs in the Recruitment module (time steps 2-30)
dplyr::rows_update(
tibble::tibble(
label = "log_devs",
time = 2:30,
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
)
),
by = c("label", "time")
) |>
# Update log_sd for log_devs in the Recruitment module
dplyr::rows_update(
tibble::tibble(
module_name = "Recruitment",
label = "log_sd",
value = 0.4
),
by = c("module_name", "label")
) |>
# Update inflection point and slope parameters in the Maturity module
dplyr::rows_update(
tibble::tibble(
module_name = "Maturity",
label = c("inflection_point", "slope"),
value = c(2.25, 3)
),
by = c("module_name", "label")
) |>
# Update log_init_naa values in the Population module
dplyr::rows_update(
tibble::tibble(
label = "log_init_naa",
age = 1:12,
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
)
),
by = c("label", "age")
)
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 tibble returned by create_default_parameters()
is
just a data frame containing specifications. Nothing has been created in
memory as of yet. To actually initialize the 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_4_model |>
initialize_fims(data = data_4_model) |>
fit_fims(optimize = FALSE)
clear()
# Run the model with optimization
fit <- parameters_4_model |>
initialize_fims(data = data_4_model) |>
fit_fims(optimize = TRUE)
## ✔ Starting optimization ...
## ℹ Restarting optimizer 3 times to improve gradient.
## ℹ Maximum gradient went from 0.00512 to 0.00032 after 3 steps.
## ✔ Finished optimization
## ✔ Finished sdreport
## ℹ FIMS model version: 0.7
## ℹ Total run time was 6.55601 seconds
## ℹ Number of parameters: fixed_effects=77, random_effects=0, and total=77
## ℹ Maximum gradient= 0.00032
## ℹ Negative log likelihood (NLL):
## • Marginal NLL= 3166.02383
## • Total NLL= 3166.02383
## ℹ Terminal SB= 1779.18153
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)
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.006489698 1
## 2 0.007156588 0.006553593 2
## 3 0.006553376 0.006572973 3
## 4 0.006501745 0.006459369 4
## 5 0.006562572 0.006211229 5
## 6 0.006467796 0.006193904 6
## 7 0.005248220 0.006068948 7
## 8 0.006568311 0.005916306 8
## 9 0.005516030 0.005303574 9
## 10 0.005614254 0.005111184 10
## 11 0.004269091 0.004836804 11
## 12 0.004842784 0.004510400 12
## 13 0.004018395 0.004146893 13
## 14 0.004741539 0.003984685 14
## 15 0.003444776 0.003705819 15
## 16 0.004062555 0.003471227 16
## 17 0.003790370 0.003369435 17
## 18 0.002892876 0.002924379 18
## 19 0.002765623 0.002643862 19
## 20 0.002349581 0.002361533 20
## 21 0.002007219 0.002130067 21
## 22 0.001883646 0.001897533 22
## 23 0.001908972 0.001950093 23
## 24 0.001596426 0.001726015 24
## 25 0.001248675 0.001575694 25
## 26 0.001459593 0.001466748 26
## 27 0.001255182 0.001404757 27
## 28 0.001314476 0.001365663 28
## 29 0.001364266 0.001292371 29
## 30 0.001343069 0.001401679 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.6443 1
## 2 461.0895 461.0768 2
## 3 747.2900 747.2702 3
## 4 996.9711 996.9800 4
## 5 767.5477 767.5611 5
## 6 1343.8609 1343.9332 6
## 7 1319.0663 1319.1035 7
## 8 2597.5051 2598.0667 8
## 9 1425.8608 1425.9444 9
## 10 1643.6742 1643.6463 10
## 11 1770.8245 1770.6432 11
## 12 1751.2628 1751.0176 12
## 13 1194.0655 1193.9453 13
## 14 1638.4036 1638.3372 14
## 15 1584.1855 1584.3797 15
## 16 1333.1820 1333.2114 16
## 17 2325.4425 2325.8491 17
## 18 1693.7841 1693.9358 18
## 19 1531.8238 1532.1477 19
## 20 1358.8162 1359.3601 20
## 21 1638.8798 1639.0245 21
## 22 1077.3713 1077.1712 22
## 23 1647.9104 1647.3001 23
## 24 1224.4634 1224.6413 24
## 25 1126.9182 1126.7708 25
## 26 956.4529 956.3227 26
## 27 902.1621 902.1436 27
## 28 1168.8217 1168.8928 28
## 29 852.3819 852.0580 29
## 30 1274.3777 1274.3621 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_4_model |>
# Update the slope value of the logistic selectivity for the survey
dplyr::mutate(
value = dplyr::if_else(
module_name == "Selectivity" &
fleet_name == "survey1" &
label == "slope",
2.5,
value
)
)
parameters_low_slope <- parameters_4_model |>
dplyr::mutate(
value = dplyr::if_else(
module_name == "Selectivity" &
fleet_name == "survey1" &
label == "slope",
1,
value
)
)
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.04636 to 0.00107 after 3 steps.
## ✔ Finished optimization
## ✔ Finished sdreport
## ℹ FIMS model version: 0.7
## ℹ Total run time was 6.36311 seconds
## ℹ Number of parameters: fixed_effects=77, random_effects=0, and total=77
## ℹ Maximum gradient= 0.00107
## ℹ Negative log likelihood (NLL):
## • Marginal NLL= 3166.02383
## • Total NLL= 3166.02383
## ℹ Terminal SB= 1779.19272
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.02164 to 0.00321 after 3 steps.
## ✔ Finished optimization
## ✔ Finished sdreport
## ℹ FIMS model version: 0.7
## ℹ Total run time was 6.28882 seconds
## ℹ Number of parameters: fixed_effects=77, random_effects=0, and total=77
## ℹ Maximum gradient= 0.00321
## ℹ Negative log likelihood (NLL):
## • Marginal NLL= 3166.02383
## • Total NLL= 3166.02383
## ℹ Terminal SB= 1779.19885
clear()
Age only
The same model can be fit to just the age data, removing the length composition configurations.
# Create default parameters, update with modified values, initialize FIMS,
# and fit the model
age_only_fit <- parameters_4_model |>
# remove rows that have module_type == LengthComp
dplyr::rows_delete(
y = tibble::tibble(module_type = "LengthComp")
) |>
initialize_fims(data = data_4_model) |>
fit_fims(optimize = TRUE)
## Matching, by = "module_type"
## ✔ Starting optimization ...
## ℹ Restarting optimizer 3 times to improve gradient.
## ℹ Maximum gradient went from 0.00513 to 0.00134 after 3 steps.
## ✔ Finished optimization
## ✔ Finished sdreport
## ℹ FIMS model version: 0.7 ℹ Total run time was 0.56636 seconds ℹ Number of
## parameters: fixed_effects=77, random_effects=0, and total=77 ℹ Maximum
## gradient= 0.00134 ℹ Negative log likelihood (NLL): • Marginal NLL= 1565.4676 •
## Total NLL= 1565.4676 ℹ Terminal SB= 1724.28719
clear()
Length
The same model can be fit to just the length data, removing the age composition configurations.
# Create default parameters, update with modified values, initialize FIMS,
# and fit the model
length_only_fit <- parameters_4_model |>
# remove rows that have module_type == AgeComp
dplyr::rows_delete(
y = tibble::tibble(module_type = "AgeComp")
) |>
initialize_fims(data = data_4_model) |>
fit_fims(optimize = TRUE)
## Matching, by = "module_type"
## ✔ Starting optimization ...
## ℹ Restarting optimizer 3 times to improve gradient.
## ℹ Maximum gradient went from 0.02233 to 0.00021 after 3 steps.
## ✔ Finished optimization
## ✔ Finished sdreport
## ℹ FIMS model version: 0.7 ℹ Total run time was 4.79269 seconds ℹ Number of
## parameters: fixed_effects=77, random_effects=0, and total=77 ℹ Maximum
## gradient= 0.00021 ℹ Negative log likelihood (NLL): • Marginal NLL= 1520.03751 •
## Total NLL= 1520.03751 ℹ Terminal SB= 1706.88667
clear()