Random Effects in Fisheries Integrated Modeling System (FIMS): Architecture, implementation, and future development

class: center, middle, inverse, title-slide

.title[
# Random Effects in Fisheries Integrated Modeling System (FIMS): Architecture, implementation, and    future development
]
.subtitle[
## FIMS CIE Review
]
.author[
### Andrea Havron, NOAA Fisheries OST <a href="mailto:andrea.havron@noaa.gov" class="email">andrea.havron@noaa.gov</a> 
]
.date[
### 2026/04/23
]

---

layout: true

.footnote[U.S. Department of Commerce | National Oceanic and Atmospheric Administration | National Marine Fisheries Service]

code.cpp{
  font-size: 20px;
}
code.r{
  font-size: 20px;
}
code.Rcpp{
  font-size: 20px;
}

</style>

---
# Outline

- Hierarchical Modeling in Stock Assessments

- Guiding Principles in Design

- C++ architecture

- User Interface

- Informing Decisions via Sparsity Analysis

- Future Development Plans

---

# Hierarchical Modeling in Stock Assessments

.center[
.large[

`$Index \sim LogNormal(\hat{\mathbf{I}}, \sigma_{I})$`

`$\hat{\mathbf{I}} = f(Growth, Maturity, Mortality, \mathbf{Recruitment}, Selectivity)$`

`$\log(\mathbf{Recruitment}) \sim Normal(\log(\hat{R}), \sigma_{R})$`
]]

---

# Fitting Hierarchical Models: Speed vs. Accuracy

.pull-left[
**Penalized likelihood**
- Fixed effects estimation (fastest)
- Uses a Constrained optimization approach

**Laplace approximation**
- Random effects estimation (moderate)
- Uses 2nd-order Taylor series approximation

**MCMC**
  - Full posterior sampling (slowest)
  - Asymptotically exact posterior
]

.pull-right[
![](slides_20260423_Random_Effects_files/figure-html/unnamed-chunk-1-1.png)
]

---
# FIMS Random Effects Requirements

- Include random effects option for fitting state-space models

- Temporal and spatial varying random effects

- Capture process variability and measurement error

- Shared parameterizations across fleets and populations

- Multivariate random effects

- Hierarchical estimation across processes, species, areas, etc.

---
# Guiding Principles in Design

- Generic and Flexible

- Extensible

- Separate the biological math from the statistics

- Dependency on TMB and access to the Laplace Approximation
---
# Choosing the default

.center[
.quote-card[
Unbiased estimation of parameter deviations (such as recruitment or selectivity) in models without random effects requires iterative bias adjustment algorithms which are inefficient and prone to error.
]
]

&emsp;&emsp;&emsp;&emsp;&emsp; ~ FIMS Requirements Documents
---
.large[**Decoupling Biological Processes from Statistical Likelihood**]
<div id="htmlwidget-1ceee326362fdfe72577" style="width:100%;height:85%;" class="grViz html-widget"></div>
<script type="application/json" data-for="htmlwidget-1ceee326362fdfe72577">{"x":{"diagram":"\ndigraph fims {\n # Global settings\n graph [rankdir=LR, bgcolor=\"transparent\", compound=true, nodesep=0.5, ranksep=0.6]\n node [shape=box, style=\"filled, rounded\", color=\"#5b8db8\", fillcolor=\"#f0f7ff\", fontname=\"Arial\", fontcolor=\"#333333\", fontsize=28, width=2.5, height=0.8]\n edge [color=\"#666666\", arrowhead=vee, penwidth=1.2, fontsize = 24]\n\n # User Layer\n subgraph cluster_user {\n label=\"User Layer\"; fontname=\"Arial-Bold\"; fontcolor=\"#5b8db8\"; style=\"dashed\"; color=\"#5b8db8\"; fontsize = 28;\n node [fillcolor=\"#e1f5fe\", color=\"#01579b\"]\n A [label=\"R User Interface\"]\n }\n\n # Central Dispatch\n subgraph cluster_logic {\n label=\"Central Dispatch\"; fontname=\"Arial-Bold\"; fontcolor=\"#c88a2a\"; style=\"dashed\"; color=\"#c88a2a\"; fontsize = 28;\n node [fillcolor=\"#fff3e0\", color=\"#e65100\"]\n B [label=\"Rcpp Interface\"]\n C [label=\"FIMS.cpp\"]\n D [label=\"Information\"]\n }\n\n # Calculation Engine\n subgraph cluster_engine {\n label=\"Calculation Engine\"; fontname=\"Arial-Bold\"; fontcolor=\"#4f9a63\"; style=\"dashed\"; color=\"#4f9a63\"; fontsize = 28;\n node [fillcolor=\"#e8f5e9\", color=\"#2e7d32\"]\n \n E [label=\"Population Dynamics\"]\n H [label=\"Catch-At-Age\"]\n F [label=\"Distributions\"]\n J [label=\"Data\"]\n }\n\n # Results\n subgraph cluster_results {\n label=\"Results\"; fontname=\"Arial-Bold\"; fontcolor=\"#b85a7a\"; style=\"dashed\"; color=\"#b85a7a\"; fontsize = 28;\n node [fillcolor=\"#fce4ec\", color=\"#880e4f\", fontsize = 24]\n G [label=\"model.hpp\"]\n I [label=\"Negative Log-Likelihood\"]\n }\n\n # Layout constraints (forcing order)\n {rank=same; B C}\n {rank=same; E H F J}\n\n # Edges\n A -> B\n A -> C\n B -> D\n D -> E\n D -> H\n\n D -> F\n D -> J \n \n J -> F [constraint=false, style=solid]\n E -> H\n \n # Grouping the outputs to model.hpp\n {H F} -> G\n C -> G\n G -> I\n \n # Feedback loops (using \"constraint=false\" so they don\"t mess up the LR flow)\n I -> C [constraint=false, style=dotted, label=\"value/gradient\", fontsize = 28]\n}\n","config":{"engine":"dot","options":null}},"evals":[],"jsHooks":[]}</script>
---
# Globally Available Variable Map

```cpp
std::unordered_map<uint32_t, fims::Vector<Type>*> variable_map;
```

uint32_t: Parameter's unique ID 
fims::Vector<Type>*: pointer to a FIMS vector

C++ map that stores pointers to all model parameters at the global scope

---
# Globally Available Variable Map

```cpp
std::unordered_map<uint32_t, fims::Vector<Type>*> variable_map;
```

.midi[
uint32_t: Parameter's unique ID 
fims::Vector<Type>*: pointer to a FIMS vector

.pull-left[
**Rcpp Interface** Populates the variable map with pointers to all model parameters, including random effects. 
**Information** Stores the variable map, which is accessible across all C++ files. 
**Distributions** Accesses the variable map to compute likelihood contributions. 
]]

---
# Globally Available Variable Map

```cpp
std::unordered_map<uint32_t, fims::Vector<Type>*> variable_map;
```

.midi[
uint32_t: Parameter's unique ID 
fims::Vector<Type>*: pointer to a FIMS vector

.pull-right[
**Prior Likelihood** 
observed_value -> parameter 
 
**Random Effects Likelihood** 
observed_value -> random effect parameter 
expected_value -> derived quantity (or vector of 0s)
 
**Data Likelihood** 
observed_value -> data vector 
expected_value -> derived quantity 
]]
---
# Connecting Random Effects to TMB

``` r
recruitment$log_devs[i]$set_estimation_type("random_effects")
```

In the Rcpp Interface, if a parameter's estimation type is set to "random_effects":
- The parameter gets added to the global `random_effects_parameters` vector.
- An R function returns the list of random effects parameters for a model.
- The random effects list is passed into FIMS.cpp through MakeADFun and is linked
   to the global `random_effects_parameters` vector.

``` r
parameters <- list(
 p = get_fixed(),
 re = get_random()
 )
 obj <- TMB::MakeADFun(
 data = list(), parameters, DLL = "FIMS",
 silent = TRUE, map = map, random = "re"
 )
```

---

# Default FIMS Model

``` r
library(FIMS)
# Prepare the package data for being used in a FIMS model
data("data_big")
data_4_model <- FIMSFrame(data_big)

# Set up the default configurations and parameters for the model
default_configurations <- 
 create_default_configurations(data = data_4_model)
default_parameters <- create_default_parameters(
 configurations = default_configurations,
 data = data_4_model
)

# run the model
fit <- default_parameters |>
 initialize_fims(data = data_4_model) |>
 fit_fims(optimize = TRUE)
```

---
# Recruitment Deviations are random effects by default

``` r
default_parameters |> 
  tidyr::unnest(cols = data) |>
  dplyr::filter(
    label == "log_devs"
  ) 
```

```
## # A tibble: 29 × 12
## model_family module_name fleet_name module_type label age length time value estimation_type distribution_type distribution
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <chr> 
## 1 catch_at_age Recruitment <NA> BevertonHolt log_devs NA NA 2 0 random_effects process Dnorm 
## 2 catch_at_age Recruitment <NA> BevertonHolt log_devs NA NA 3 0 random_effects process Dnorm 
## 3 catch_at_age Recruitment <NA> BevertonHolt log_devs NA NA 4 0 random_effects process Dnorm 
## 4 catch_at_age Recruitment <NA> BevertonHolt log_devs NA NA 5 0 random_effects process Dnorm 
## 5 catch_at_age Recruitment <NA> BevertonHolt log_devs NA NA 6 0 random_effects process Dnorm 
## 6 catch_at_age Recruitment <NA> BevertonHolt log_devs NA NA 7 0 random_effects process Dnorm 
## 7 catch_at_age Recruitment <NA> BevertonHolt log_devs NA NA 8 0 random_effects process Dnorm 
## 8 catch_at_age Recruitment <NA> BevertonHolt log_devs NA NA 9 0 random_effects process Dnorm 
## 9 catch_at_age Recruitment <NA> BevertonHolt log_devs NA NA 10 0 random_effects process Dnorm 
## 10 catch_at_age Recruitment <NA> BevertonHolt log_devs NA NA 11 0 random_effects process Dnorm 
## # ℹ 19 more rows
```

---
# Change to fixed effect
Estimate with constrained optimization (penalized likelihood)

``` r
# set estimation_type of log_devs to fixed_effects
default_parameters <- default_parameters |>
 tidyr::unnest(cols = data) |>
 dplyr::mutate(
 estimation_type = ifelse(label == "log_devs", 
 "fixed_effects", estimation_type)
 )
# fix Recruitment Distribution log_sd to "constant"
default_parameters <- default_parameters |>
 dplyr::mutate(
 estimation_type = ifelse(label == "log_sd" & 
 module_name == "Recruitment", 
 "constant", estimation_type)
 )
```
---
# Upcoming User Interface Enhancements

New helper functions to set random effects for any process (GitHub Issue [#1235](https://github.com/NOAA-FIMS/FIMS/issues/1235))

``` r
add_process(
  module_name = “selectivity”,
  fleet_name = “fleet1”,
  specification_type = “semi_parametric”,
  strucutre = “AR1”,
  base_function = “logistic”,
  logit_rho = 0)

add_priors(
    fleet = c(fleet1, fleet2),
    module = selectivity,
    slope ~ normal(1.5, 10),
    inflection_point ~ normal(2, 10)
)
```

---
# Advanced User Implementation

``` r
# Create a recruitment module
recruitment <- new(BevertonHoltRecruitment)

# Set the log_devs parameter to be estimated as random effects
recruitment$log_devs$resize(n_years - 1)
recruitment$log_devs$set_all_random(TRUE)

# Set up a recruitment distribution module
recruitment_distribution <- methods::new(DnormDistribution)
recruitment_distribution$expected_values$resize(n_years - 1)
recruitment_distribution$log_sd[1]$set_estimation_type("fixed_effects")

# Link the observed values of the recruitment_distribution
# to the log_devs parameter vector
recruitment_distribution$set_distribution_links("process", 
  recruitment$log_devs$get_id())
```
🚧 : The `set_distribution_links` function is being refactored in GitHub Issue [#1194](https://github.com/NOAA-FIMS/FIMS/issues/1194)

---
# Sparse parameterization of recruitment

$$
`\begin{align}
log\_r &\sim Normal(\widehat{log\_r}, \sigma_{r})\\
\widehat{log\_r} &= log(BevertonHolt(SSB, parameters))
\end{align}`
$$

``` r
# Create a recruitment module
recruitment <- new(BevertonHoltRecruitment)

# Set the log_devs parameter to be estimated as random effects
recruitment$log_r$resize(n_years - 1)
recruitment$log_r$set_all_random(TRUE)

# Set up a recruitment distribution module
recruitment_distribution <- methods::new(DnormDistribution)
recruitment_distribution$log_sd[1]$set_estimation_type("fixed_effects")

# Link the observed values of the recruitment_distribution
# to the log_r parameter vector and log_expected_recruitment derived quantity
recruitment_distribution$set_distribution_links("process", 
  c(recruitment$log_r$get_id(), recruitment$log_expected_recruitment$get_id()))
```
---
# Hessian Inversion with the Cholesky Factorization

To avoid the computational burden of a full matrix inversion, we use the Cholesky Factorization

$$
H = LL^T\mathrm{~where~} L \mathrm{~is~a~lower~triangular~matrix}
$$

.pull-left[.large[**Penalized Likelihood**]
 
.left[
- Uses the Cholesky factor, `$L$`, to calculate the Variance-Covariance matrix **once** at the MLE for uncertainty estimation

]
]

.pull-right[.large[**Laplace Approximation**]
 
.left[
- Uses the Cholesky factor, `$L$`, to evaluate the Marginal Likelihood during **every iteration** of the optimizer
]]
 
**This factorization is efficient but its performance depends entirely on matrix structure**
---

# The Scalability Bottleneck

.large[The ] `$\Large{O(n^3)}$` .large[ wall]

.large[
- Dense Cholesky factorization scales at `$O(n^3)$`
 
- In the Laplace Approximation this cost is paid at **every iteration**
 
- To move beyond moderate-sized models, we must exploit **Sparsity**]

---
# Sparsity and TMB
- Laplace Approximation requires the **Cholesky factorization** of the Hessian at **every iteration**
 
- Template Model Builder (TMB) makes this fast by taking advantage of **sparsity**
 
- After building the **static** computational graph, TMB **automatically detects** sparseness and **optimizes** the tape by removing operations involving constants of `$0$`
 
- When evaluating the Laplace Approximation, TMB uses the **sparse Cholesky factorization**, which eliminates all `$0$` calculations
 
`$O(n^3)$`: Dense Hessian 
`$O(n^{3/2})$`: Sparse Spatial Hessian (GMRF / 2D Grid) 
`$O(n)$`: Sparse Time Series Hessian (AR1, RW / 1D)
---
# Parameterization affects sparsity of Hessian

.mylarge[
**AR1 model:** `$x_t = \phi x_{t-1} + w_t$`]

.pull-left[
**Process Parameterization**

`$x_t \sim N(\phi x_{t-1}, \sigma^{2}_{x})$`

Sparse Hessian, Cholesky factorization: `$O(n)$`
]

.pull-right[
** Deviations Parameterization**

.narrowtopbottommargin[
`$x_{t} = \phi x_{t-1} + \sigma_x w_{t}$` 
`$w_{t} \sim N(0,1)$`]

.normal[
Dense Hessian, Cholesky factorization: `$O(n^3)$`]
]

`$w_t$` is white noise and `$\phi \neq 0$` for an order-_1_ process

---
# Benchmarking AR1
.pull-left[
 Speed 
<img src="static/ar1_benchmark.png" alt="Benchmarking the speed of AR1 Parameterizations" width="100%" />
]

.pull-right[
 Memory 
<img src="static/ar1_memory.png" alt="Benchmarking the memory of AR1 Parameterizations" width="100%" />
]
---
# Applications in Stock Assessments

.pull-left[
**Case Study: Catch-at-age Model**

- Models fit using a modified **babySAM** in RTMB
 
- Recruitment, Numbers at Age, and Fishing Mortality are AR1 processes
 
- **Focus:** Comparing parameterizations in Recruitment
]

.pull-right[
**Parameterizing Recruitment**

**Process:**
`\begin{aligned}
log(\widehat{N_{y,1}}) &= log(N_{y-1,1})\\
log(N_{y,1}) &\sim Normal(\bar{r} + \phi(log(\widehat{N_{y,1}}) - \bar{r}), \sigma^{2}_{r})
\end{aligned}`

**Deviations:**
.narrowtopmargin[
`\begin{aligned}
z_{1} &\sim Normal\big(0, \sqrt{\sigma_{r}^{2}/(1 - \phi^2)}\big)\\
z_{2:n} &\sim Normal(0, \sigma_{r})\\
dev_{1} &= z_{1}\\
dev_{y} &= \phi dev_{y-1} + z_{y}\\
log(N_{y,1}) &= \bar{r} + dev_{y}
\end{aligned}`
]]

---
# Speed and Memory Test Results

- Relative results from 100 interactions using `bench::mark`
- Model fit `$n=45$` years of recruitment

|Intercept|Model       |Relative time difference |Relative memory allocation  |
|---------|------------|-------------------------|----------------------------|
|Y        |**process** |**1**                    |**1**                       |
|N        |**process** |**1**                    |**1**                       |
|Y        |deviations  |1.4                      |1.4                         |
|N        |deviations  |1.6                      |1.4                         |

At these sample sizes, the ~**1.5x speedup** reflects the reduction in the recruitment block's computational cost. The true `$O(n^3)$` "wall" becomes more visible as model scales up in complexity

---
# Moving average (MA) models

**MA1 model:** `$x_t = w_t + \theta w_{t-1}$`
.pull-left[
**Process Parameterization**

`$w_{t} = x_{t} - \theta w_{t-1}$` 
`$w_t \sim N(0, \sigma_{x})$` 
`$y_{t} \sim N(x_{t}, \sigma_{y})$` 
where `$w$` and `$x$` are treated as random

]
.pull-right[
** Deviations Parameterization**

`$w_t \sim N(0, \sigma_{w})$` 
`$x_{t} = w_{t} + \theta w_{t-1}$` 
`$y_{t} \sim N(x_{t}, \sigma_{y})$` 
where `$w$` is treated as random

]
---
# Next Generation Models
 
.large[
- **Recommendation**: 
 - **AR1**: process approach
 - **MA1**: deviations approach
 
- No universal best parameterization but **optimal sparsity** depends on model structure
 
- For next gen TMB models, **good practices in efficiency and inference** should be evaluated through the lens of the Laplace Approximation
]
---
# Future FIMS RE Development Goals
 
.large[
- Add a generic **ARn** distribution (Issue [#582](https://github.com/NOAA-FIMS/FIMS/issues/582))
 
- Add multivariate functionality 
 - **GMRF** distribtion (Issue [#582](https://github.com/NOAA-FIMS/FIMS/issues/582))
 - Ability to set **multivariate pointers**
 
- Fully develop **non-parametric** and **semi-parametric** functionality for all time-varying parameters 
 - Age-based selectivity Issue [#1206](https://github.com/NOAA-FIMS/FIMS/issues/1206)
 
- Refactor **Numbers At Age** with random effect capability
 
- Add OSA residuals (Issue [#546](https://github.com/NOAA-FIMS/FIMS/issues/546))
]
---
# Summary
 
.large[
- FIMS Random Effects are designed to be **generic**, **flexible**, and **extensible**
 
- FIMS keep the biological math **separate** from the statistical likelihood 
 
- The global **variable map** allows any **parameter** or **derived quantity** to be linked to random effect distributions
 
- The current **default model** fits recruitment deviations as **random effects** with a Normal iid distribution
]

---
# &#x1F44B; Thank you!

&#x1F4E9; andrea.havron@noaa.gov